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ABSTRACT 

Two approaches to solving problems involving thin 
shells based on the standard methods of structural analysis 
are discussed. In the stiffness method, the governing shell 
equations are expanded into a Fourier series and reduced to 
a set of eight first order differential equations. A forward 
numerical integration technique is used to form the 
stiffness matrix and the particular solutions. In the 
flexibility method, the governing shell equations are 
Simplified by limiting the analysis to axisymmetric shells 
of constant thickness. Closed form solutions are obtained 
for the flexibility coefficients for specific shell 
geometries. Particular solutions are approximated by the 
appropriate membrane solution. 

A computer program was developed to perform the 
flexibility analysis based on the approach presented. The 
results are compared with the results from a program 
developed by Shazly (5) based on the stiffness method. The 


solutions from the two programs show excellent agreement. 
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NOMENCLATURE 
radius of curvature of a sphere 
matrix coefficients which are a function 
of the geometric and material properties 
of the shell 
Boolean Connectivity Matrix 
load vector coefficients 
constants of integration vector 
constant parameter which is a function 
of the rigidities of the shell and the 
principal shell curvature 
flexural rigidity 
segment deformation vector 
modulus of elasticity 


fixed end forces 


primary and secondary force vectors 
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segment flexibility matrix 
Structure flexibility matrix 

shell thickness 

horizontal force, positive in the 
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revolution 

£Ficthtiondlahonazontalt forcej duefto 
a vertical edge load at the top of 


a cone or sphere 


homogeneous solution vector used in the 
stiffness analysis 


transfer matrix arising from integrating 
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linear differential operator defined 
in Eqn. A.11 


in-plane bending moments 


twisting moments 
harmonic number 


normal in-plane forces 


in-plane shear forces 


intensity of the load components in the 
directions s,¢,6,z respectively 


particular solution vector used in the 
stiffness analysis 


structure particular solution used in 
the flexibility analysis 


vector arising from the integration 
OnewBe} 


coefficients of {0,} 
transverse shear forces 


radius of the parallel circle for the 
cylindrical segment 


radius of a parallel circle 
radius of curvature of a meridian 
length of the normal between any point 
on the midsurface and the axis of 
revolution 

curvature of a parallel circle 

first principal curvature=1/r, 

second principal curvature=1/r2 


total vertical load acting on a segment 
due to the applied loads 


coordinate which measures the distance 
along the shell meridian 
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effective transverse shear force 
effective tangential shear force 


matrix relating the constants of 
integration to the redundant vector 

for the segment; a function of the shell 
geometry 


matrix relating the constants of 
integration to the particular solution 
displacement vector; a function of the 
geometric and material properties of 
the shell 


displacement component in the 
circumferential direction 


change of variable in terms of Qg, and 
r2 used to form the homogeneous 
solution in the flexibility analysis 


displacement component in the 
meridional direction 


change of variable in terms of r,,Vv,w, 
used to form the homogeneous solution 
in the flexibility analysis 


displacement component in the radial 
direction 


vector of the eight dependent variables 
in the stiffness analysis 


coordinate which measures the distance 
in the direction normal to the 
midsurface toward the axis of revolution 


angle between the outer edge of the 
sphere and the axis of revolution, or 
the semi-vertex angle for a cone 


angle between the inner edge of the 
sphere and the axis of revolution 


coefficient of thermal expansion 
parameter which is a function of p, 
r, and h in the flexibility analysis, 
or, the meridional rotation in the 
stiffness analysis 

specific weight of the shell or the 
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liquid weight density 
shear strain 


particular and homogeneous deformations 
respectively 


horizontal displacement of a shell 
meridional rotation of a shell 
change of variable used to form the 
homogeneous solution for the cone, 


defined in Eqn. A.19 


coordinate which measures the angle in 
the circumferential direction 


dimensionless parameter which is a 
function of a/h and’ yp? for ‘the 
Sphere, or the parameter in terms of 
h ane the semi-vertex angle for the 
cone 


Poisson's ratio 

dimensionless input parameter for the 
evaluation of the Kelvin® functions for 
the cone 

meridional and circumferential stresses 
Shear stress 

coordinate which measures the angle 


between any point on the midsurface and 
the axis of revolution 
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1. INTRODUCTION 


1.1 Introductory Remarks 

A shell of revolution is a surface generated by 
rotating a plane curve about an axis lying in the same 
PLranetashebiscotyuevolutnon fornmivpartrof such structures. as 
pressure vessels, storage tanks, silos, nuclear containment 
Structures, and cooling towers. Apart from their attractive 
appearance, the widespread use of such shells as structural 
elements is attributed to their efficiency in resisting 
load. This leads to thinner sections and reduced material 
costs. 

The general theory of shells of revolution, originally 
developed by Flugge, applies to any type of meridian 
geometry with either constant or variable thickness, and 
Subjected to any type of loading. However, for many 
practical applications the shell segments are of constant 
thickness and the loads are axisymmetric. Tegan are of 
such shells can be simplified by separating the solution 
into two partssefinsuly, Gtheepanrtieulanysolution 
approximated by the membrane stresses due to the applied 
loads: and secondly, the bending stresses due to the edge 
effects. Moreover, if the shell segments are sufficiently 
long such that there isevirtually no interaction between the 
edges of a shell segment, the computations can be simplified 
even further. This method is analogous to the method of 


consistent deformation in elastic frame analysis. And since 
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accounting for the boundary effects involves evaluating the 
flexibility coefficients, this method of analysis will be 


referred tovasthe flexibility method: 


1.2 The Objectives of the Study 
The objectives of this study are: 

1. To review the solutions to the general theory of shells 
of revolution; 

2. To obtain solutions for the membrane stresses and 
flexibility influence coefficients in closed form for 
cylindrical, spherical, and conical segments under 
various axisymmetric loadings. 

3. To incorporate these solutions into the computer program 
FLEXSHELL. 

4. To evaluate the limitations of an approximation used in 
obtaining the solution for spherical segments known as 


Geckeler's assumption. 


1.3 Structure of Thesis 

The thirteen basic differential equations of shells of 
revolution are formulated in detail in Chapter 2. Chapter 3 
presents the two solution techniques to solve these 
governing shell equations based on standard methods of 
structural analysis. The formulation of program FLEXSHELL 
based on the flexibility approach is presented in Chapter 4. 
An evaluation of the accuracy of the closed form solutions 


used in this approach is presented in Chapter 5. Finally, 
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Chapter 6 consists of a brief summary and conclusions of the 
study. Detailed derivations, the program listing, sample 
input and output files, and the user's manual for program 


FLEXSHELL are found in the Appendices. 
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2. THEORY OF SHELLS OF REVOLUTION 


2.1 Shell Geometry 

As shown in Fig. 2.1, a shell is geometrically defined 
by its midsurface which bisects the shell thickness, h. A 
surface of revolution is generated by the rotation of a 
plane curve about an axis in its plane. This generating 
curve iS called a meridian. Another term frequently used is 
the parallel circle, which is the intersection of the 
surface with a plane perpendicular to the axis of 
revolution. To specify an arbitrary point on the midsurface, 
two coordinates need be specified: 6, the angular distance 
of the point from the datum meridian, and ¢, the angle 
between a normal to the shell and its axis of revolution. To 
measure the distance along a normal to the midsurface, a 
third coordinate z, may be specified. The radii of curvature 


of a shell of revolution are: 


Pe = radius of the paratled”circle; 
re.=tradius of curvature of 4 meridian: 
r2 = length of the normal between any point on 


the miadsurface and the axis of revolution. 


The following relations can be derived from Fig. 2.1. 
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The internal stress resultants in Fig. 2.2, is 
determined by integrating the internal stresses through the 


shell thickness as follows 
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The Fundamental Assumptions 


The fundamental equations of the general theory of 


shells of revolution first presented by Flugge (1) are based 


on the following set of assumptions: 


lee 


Now, 


Thin shell - the shell thickness is small in comparison 
to the other dimensions of the shell. Thus, the stresses 
on the z-face, and the twisting moments about the z-axis 
may be neglected. 

Small deflection theory applies. The displacements of 
the shell due to the applied loads are sufficiently 
small that the equilibrium equations developed from the 
initial shell geometry do not change. 

Material is linearly elastic, i.e., Hooke's law applies. 
Plane sections remain plane after bending. i.e., the 
normals to the middle surface before bending remain 
normal after bending. 

Deformations of the shell due to radial shears can be 
neglected. 


based on these set of assumptions and the shell 


geometry, the general theory of shells of revolution may be 


formulated by: 


Pe 


Determining the equilibrium of forces acting on the 
Gutferential element shown in Fig. 222° (six ‘equations 
with ten unknowns) 

Establishing the strain-displacement relationships; (six 
equations with six unknowns) 


Establishing the stress-strain relationships from 
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Hooke's Law; (three equations with six unknowns) 

4. Transforming the stress-resultant equations into the 
force-displacement equations; (six equations with three 
unknowns) | 

5. Obtaining a complete formulation by combining the 
force-displacement equations with the equilibrium 


equations. (thirteen equations with thirteen unknowns) 


2.3 Equations of Equilibrium 

Conciderjine, dipierential element shown Jngrig. 2.2. 
From the summation of forces in each of the coordinate 
directions and moments about each of the coordinate axes, on 


6, and z, the six equations of equilibrium are: 
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No cue Newaseieridi ones andecircimiterential shear forces; 
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Qo, Qe = transverse shear forces; 
Mo, Me = meridional and circumferential moments, 
respectively; 


Moo, Moo = meridional and dieeane zenei at twisting 
moments, respectively. 

Note that all forces and moments are expressed in units 
of force per unit length. The sign convention used is as 
shown in Fig. 2.2, where Ng and Ne are positive for tension 
along the meridian and ay Menneer Once vesnect velit My and 
Me are positive when the Bice: shell surface is in 


compression. 


2.4 Force-Displacement Equations 

The deformation of a shell element consists of the 
change in length of the shell edges, r,d¢ and rod6, and of 
the change of the angle between these edges. In reference to 
Fig. 2.3, the midsurface strain-displacement relationships 


for a shell element are: 


Meridional strain, eg = 1 (v’- w) 2.4(a) 
Tr; 
Hoop strain, €o = 1 (u' + veos¢@ - wsin®¢) 2.4(b) 
To 
Shear strain, Noe = Ne) = Ue UECOS¢ 2.4(c) 
Yo a4 To ; 


where 

u = midsurface displacement component in the circumferential 
direction, positive in the direction of increasing 6. 

v = midsurface displacement component in the meridional 


direction, positive in the direction of increasing ¢. 
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w = midsurface displacement component in the radial 
direction, positive in the direction away from the centre of 
curvature. 

Consider a point i at a distance z to the midsurface, 
een Ol a Soba + Sy and (ro) | =) ree tez,. FromaEqn.. 2..1(a); 


the strains at point 1 are: 
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Hooke's law forms the basis for the formulation of the 


stress-strain equations. 
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where E is the modulus of elasticity and v 1s Poisson's 
ratio. Combining the strain-displacement relationships 
(Eqns. 2.5 and 2.6) and’ substituting these into the 
stress-strain equations (Eqns. 2.7), and finally, 


substituting these into the stress-resultant equations 
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(Eqns. 2.2) and integrating through the shell thickness, the 


force-displacement relationships are as follows: 
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Where the extensional rigidity K and the flexural rigidity 
D, are defined respectively as 


K = Eh 
(t-v?) 


Bh? 
(1-p?) 


There are now fourteen equations (Eqns. 2.3 and 2.8) 


0 
" 


with thirteen unknowns, Ng, No, Nog, Noo, Mo, Mo, Moa, Moo, 
Qo, Qe, u, Vv, Ww. Note that there is one equation too many. 
Since both sides of Eqn. 2.3(f) are small differences 
between small quantities which are almost equal, this 
equation may be discarded. Thus, there is now a balance of 
unknowns and equations. The classical method of solution 
would be to reduce these differential equations into a 
Single eighth order equation in terms of one variable. This 
procedure tends to be too complicated and cumbersome to 
solve. Therefore, alternative solutions to these equations 
based on the standard methods of structural analysis will be 


presented in the following chapter. 
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AXIS OF REVOLUTION 
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Figure 2.1 GEOMETRY OF SHELL 
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Figure 2.2 FORCES ACTING ON 
SHELL MIDSURFACE 
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Figure 2.3 SHELL SEGMENTS BEFORE 
& AFTER DEFORMATION 
(a) meridian 
(b) parallel circle 
(c) angle change 
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3. METHOD OF ANALYSIS 

Structures that geometrically consist of several 
segments of shells of revolution can be analyzed by either 
of the two standard methods of structural analysis, namely: 
the stiffness method and the flexibility method. With the 
stiffness method, the stiffness matrix relating the forces 
and deformations at the edge of each shell segment are 
computed using the procedures outlined in Section 3.1. These 
element stiffness matrices are then superposed to form the 
Structural stiffness matrix from which the segment edge 
deformations are computed. With the flexibility approach 
(Section7302)yathe flexrbility influence*ecoefficients for 
each shell segment are obtained. Equations of geometric 
compatibility at the segment boundaries are written to 
obtain the forces at segment junctions. 

In this study, the use of the flexibility approach is 
explained in detail for structures with axisymmetric 
loading. The membrane solutions are used for the particular 


SOLUELON. 


3.1 The Stiffness Approach 

Program SASHELL analyzes a segmented shell structure 
based on the stiffness method. To establish the stiffness 
matrix and fixed end forces vector the basic shell equations 
must be solved numerically. But in order to do this, the 
shell equations must be reduced to a set of eight first 


order differential equations, corresponding to the eight 
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natural boundary conditions of the shell segment by: 

1. Expanding the equations using a Fourier series to 
eliminate the necessity of forming the equilibrium 
equations in the circumferential direction. 

2. Introducing the auxilliary equations to eliminate the 
in-plane shear force in the circumferential direction 
and the meridional twisting moment. 

3. Performing matrix operations to eliminate the forces in 
the circumferential direction. 

Introducing the s coordinate, which measures the 
distance along the shell meridian, the five independent 


equations of equilibrium (Eqns. 2.3) become 
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The force-displacement equations (Eqns. 2.8) in terms of the 


S-coordinate are 
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3.1.1 Fourier Series 
For any variable, say F(x,y) being arbitrary functions 


of x and y, may be represented in the form 


© © 
Fe= 2 Fa (x)cos(ny) + = Fa (x)sin(ny) Sas 
where n is the harmonic number and variable F, is now a 
PUNCE TON Of ox Only. “Similarly, the load components p,, pe, 
icp?) anal Orcas eN |, No, Neo Now, Mon MoyncMe. pp Meo, Os, 
and Qe, and displacement components u, v, and w, may be 
expressed aS a Fourier series, where the variable components 
become a function of s only. The first and second series in 
each expression represent the portions of the variables 
which are respectively symmetric and anti-symmetric with 


respect to the meridian passing through the line 6 = 0. 
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subscript n (Eqn. 3.4) can be substituted into the basic 
Shell equations (Eqns. 3.1 and 3.2), because the sequences 
sin(n@) and cos(né@) are linearly independent. 
Differentiations with respect to @ can be performed and the 
terms grouped according to the common factors, cos(n@) and 
Sin(n@). Since the coefficient of each of these factors must 
be zero, each factor produces a Separate equation. For 
example, for any n, the cosine terms in Eqn. 3.1(a) become 


roN?,cos(né@) +. cos@N;,cos(né@) + nN,,),cos(né) 


= GOSGN,nCOS( ng) ->9reO7,cos (ne) + Rop.,cosing) «= Owesrs 
tT 4 
which, upon factoring out the common term, yields 
BoNicn + COSO@N,, + NNo;n — COSONG,: — oes CE Se Ol at Ge 0 3.6 
es 
Similarly, for the sine terms, Eqn. 3.1(a) become 
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Let Ro, Ri, Rz be defined as shell curvature, i.e. 


Ro = a 
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R, = at 
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Thus, for the nth set of equations, the five independent 
equilibrium equations derived from Eqns. 3.1 become 
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Méne bERgcoséMpnst nRéMoed scRecdsoMs é 6-104, coal 3.8(d) 
Mpawin i+ TRocoSséM is SUF uNRoMe dat EREGOSéMs Ke meOedn= 0 3.8(e) 
and the eight force-displacement equations obtained from 
Eqn. 3.2 become 
Danes DOR CRis8>) epudee ane |D(R,-Rz)18,° + (K(R,t+eR2) + 
DRPCRG RA) wie LeKRocos¢é tarDR? GRicRa) ro: 2ivalt ibKet 
DRaGR, SE 2h) vee lt ( benkRalus 3.94) 
Non = [DRo(Ri-R2)cos¢]B, + [K(R2+yR,) + D(R,-R2) (R3n?-R3) Iw, 
t+elKRecosé - DRER-cosé(Ri-Ri) Ivact [vuK]vae° c+ [nkRolu; 
not b) 
Negne= 05561-») 1£[nDR6 (RacRe) 1B. t [nDRicose(RisR2)]wa ct 
[nKRo + nDRoR,(Ri-Rz2)]v, - [KRocos¢ - DRocos¢(R,-R2)?Ju, 
+hhKet IDCR{oR> kaduent 32 9(Cu 
Nosn = 0.5(1-v){t{nDRo(Ri-R2) 18, + [nDRécos(&v>(Ri-R2) Jw, F 
EnKRoat inDRoRTGR Rallye =alkRecosdé)uy teobKlug? 3n9(d) 
MaresiiDR upeirivDRecesddfianslDisas £ebDRabRiszReh > 
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defined in the following section. Note that there are two 
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sets of equations, grouped according to the cosine and sine 
terms. The final solution is obtained by solving each set 


Separately and superimposing the two solutions. 


3.1.2 Auxilliary Equations 

The quantities in the natural boundary conditions on 
the edges of a shell segment are the four displacement 
components, the rotation of the meridian (8), the radial 
displacement (w), the meridional displacement (v), and the 
circumferential displacement (u), and the four corresponding 
forces, the meridional moment (M,), the effective transverse 
Shear force (S,), the normal in-plane meridional force (N,), 
and the effective tangential shear force (T,). Three of 
these variables 8, T,, and S, do not appear in the basic 
Shell equations. They may be introduced by setting up the 
so-called auxilliary equations which express these variables 
in terms of in-plane shear forces in the circumferential 
direction and the meridional twisting moment. 

Consider the side view of the top edge of the shell 
element shown on Fig. 3.1 with two adjacent elements of 
length ds = rod@. (Note that ds = r2,d@ for very small ds) 
Thesmoments acting on the infinitesimal element!ds can-be 
replaced by a set of statically equivalent forces F, and F, 
5), SuehSthat 
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Lpansverse torce O;, and the in-plane shear force, N,¢, 
respectively, yields and expression for the Kirchhoff 
shears), Sand 7; 

5, = OF Pr RoMi 

Lia oul Tito <6 
Expanding these into a Fourier series yield 

Sen = Oe, 2 RGM. oC, Sheek 6 

Ton = Ne@y ~ R2Mion Seal 
Using the geometrical relations in Eqns. 2.1 and 3.1(f£), the 
derivatives of these forces with respect to the coordinate s 
may be written as 

Su. 2 One: etn keMs. beth co  COSOM aa. Sri 

Tea P= N ee 1 io M oo hens BR GR y=) COLOM oc See es: 
Also, by superimposing Figs. 3.2(a) and (b), the angle by 
which an element of the meridian rotates during deformation 
may be expressed in terms of the displacement components as 
follows, 


baa Ween Ra 3.14 


3.1.3 Reduction of the Shell Equations 

Rewriting Ban. 3.10 to form) an expression. for O.,, “and 
substituting this into Eqns: 3.8fa) and, (d)@respectively, 
yields 
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Rewriting Eqns. 3.11, 3.13, and 3.8(e) to form 
GxDressions fOr Noo, , Nios + and Qon, Yields 

Neon = Ten + R2Mson 

Neon o= iene teRoMeon + Ro(h eRe) COLOM. >. 

Qen = Mson + RoCOSOMson * NRoMon + RoCOSbMosn 
Substituting the above expressions into Eqn. 3.8(b), and 
using the relation, Rossing = R2, yields 
Ten = -Rocosd(Ri-Rz2)M.son — RoCOSOT;n + NRoNon — RoCOSONosn 

+ NRoR2Men + RoR2COSOMosn - Don Sha y74 

Panay FEWrDVEinoungn. sil2icto form an expression for 
Qsn , and substituting this, in addition to the expressions 
for Qen and Q,, derived earlier, into Eqn. 3.8(a), gives 
Sen = ~R2Non — RiN,n + N?R3Mon F NRZCOSH(Mosn*M.on) 

RoCOSOS.n +. Den i hea | 
Eqns. 3.15 to 3.18. may be written symbolically as 
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Santa = Bein none Mean Non Neen bon: 


OUrParLNamaek Ke LOrM- 


(Eo =aiBs ceo + {B3} 3.19 
Fe 
where 
<PeSree<Mecmricge Notre a> 
<P oa Mis ae oe ae CNG et een 
<Fo> = <Mon Mosn Mson Non Nosn? 
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the load vector. These matrices are defined in Table 3.1. 
The plus and minus signs relate to the two sets of 
equations, grouped according to the cosine and sine terms in 
the Fourier series expansion. 

To form expressions for the displacement variables, 
manipulate the force-displacement equations as follows 

het 

CA, = K + DR,(R,-R2) 3.20(a) 
CAV=OK 1 DR2 (RR: } Bi 2 OGD) 

Multiply Eqn. 3.9(a) iby (Ri-R2); 

Mie Ra) = [DRACR Re) ere 1p, — [DIR aR) 1p, + 
[KURGtPR2) © DR7(RG-R2)1(Ry-R3)w, + [vDRocos¢ - 
DRiGRaakz) LOR VoRS Van tel CAgAR{oRo) lV Gt 
[yKnR,(R,-R-) Ju, 

andemultiply Eqn. 3.9(e)) by CA,/D, 

CAG Mea D4= (Roc, ap RoCOSONCAG Ew eChdpn } ULRa Rare as 
PHORSICH Wem CA Rot eT CAG (RG—Ro) lve ot 
fynRoR>CA, ]u,, 

Subtracting» the, first. from.the second expression, and 

Simplifying by means of Eqns.- 3.20 yields, 

B,° = {-CA,M,,/D +, (RizR2)Nennt [Ror n1CA2,.- BRocGoséCa y)8, - 
[CA,vn?R3Z + vKR2(Ri-R2)]w, - [vKRocosé(Ri-R2) + 
Bir, CAs IVs tel yDRoRiCA2 luni ZCA> Susie | 

Similarity. subtracting Bon 3.9(4) trom the woroduct of 


(R,-R,) and Eqn. 3.9(e), and simplify the expression using 


Eqns. 3.20 -yields 
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[yDn?R3(R,-Rz2) + R,CAz + vKR2]w, - [vKRocos¢]v, F 
[ynRoCA2z]u,}/CA2 3222 
Rewriting Eqn. 3.14 yields 
Wn’ = VaRi — Ba S23 
Binally, (substituting Eqn. 3.9(g) into 3.11, and rewriting 


the equation to form an expression for N,on,, then 


substituting this into Eqn. 3.9(c), and simplifying, 


ae {2T,,/(1-v)_* [DnRo(R,-3R2)]£, + [DnRgcos¢(R,-3R2) ]w, 
+ [nCA,Ro - DnRoR,R2]v, + [Rocos¢CA;]u,}/CA;3 3.24 
where 
CA; = K + D(R?-3R,R2+3R3) 6s) 


Eqns. 3.21 to 3.25 may be written symbolically as 


Be. = Fea Bag Winks Vir ne ene) 
Wine = Boe ven vn 

° 
Vin = Bee (UBing Wap wn Unie nie ca) 
nae = Por (BaeWnavVarlne cee 


Ome in) machi xXe-LOrn, 


{Des LAS eae aco 
Be 


where {D} and {D°} consists of displacements £,, Wn, Un, and 
v,, and their derivatives with respect to the coordinate s, 
respectively; [A, Az] is a function of the geometric and 
material properties of the shell, defined in Table 3.2. 
Again, the plus-minus Signs relate to the set of equations, 
grouped according to the sine and cosine terms in the 
Fourier series. 

The vector <F,> which appears in Eqn. 3.19, is formed 


by writing the force-displacement equations in the following 
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SEdeGm 3.9 hye Ch) ,.1(0),2 (bd pa). Inomabr xe corm, 
feet sabe solid jr) 
D 


where {D} and {D°} are defined as before. The cofficients of 
[B, Bs] is a function of the geometric and material 
properties of the shell, defined in Tables 3.3. 
iSubstitutings Eqns, 36/27, into. at9 tyields 
Pee oh =e [BadPade sed Bold [By lipthy +t Bahibhiat ibs} 
Simplifying, 
JERE he = MAD toll ALE ah Se UBiss Sig 225) 


where 


[A3] [Bo2J{B,](A,] + (B2](Bs] 

(Aue =ore. sicre(es IB, J LAG 
Combining Eqns. 3.26 and 3.28 to form a single matrix 
equation yields 

Dix A, A2]|D 0 

eee A3 Ay F. B3 | 
Matrix equation 3.29 relates, at any point, the eight 
fundamental dependent variables, that appear in the natural 


boundary conditions of shells of revolution, and their 


derivatives with respect to the independent variable s. 


3.1.4 Solution of the Governing System of Equations ~ 

To establish the stiffness matrix, the eight first 
order differential equations expanded into a Fourier series, 
represented by matrix Eqn. 3.29, must be solved numerically. 


Ine general; Sant 3.29 can be written as 
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yneteebilant fl yatesraeig 3280 
where iy,} and {y,°} are vectors of eight dependent 
variables, four displacement components and four 
corresponding forces, and their derivatives, respectively. 
[A,] is the coefficient matrix relating the variables and 
their derivatives consisting only of functions of the 
material properties and geometry of the shell. {B,} is a 
function of the applied loads. 

The general solution of Eqn. 3.30 consists of two 
parts: the homogeneous solution and the particular solution. 
From Eqn. 3.30, the form of the homogeneous part is 

Uh, ©} = ARC, § 3,00 
and the particular solution is 

UPAGNOEIA OR PLS re LEE BIS 2 

Now, consider the solution of Eqn. 3.31 for a segment 
in the region j 2 s 2 i. Let the eight arbitrary constants 
of integration be the eight boundary conditions at edge i, 
and denote these values by {C}, then 

{h,} = ic} 8.88 
SubStitewed NGwhon'9:3). sae uimitowls .3a), Pttiomets ey 15 

thie. eal Ad ACh 3.34 
Integrating this numerically, as an initial boundary value 
problem, allows the value of h, at any point in the region 
to be determined as 

hot eet, JAH ahs She 
where [H,] represents the matrix arising from the 
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arbitrary constants of integration. For Egn. 3.35 to reduce 
LO Eqn. 3.33 when s =a, [Hi] must be the identity matrix, 
eel, 

(H,;] = {1] 3'.36 
Eqn. 3.36 may be considered to be a 'boundary condition' on 
the numerical integration of [H,]. 

Now, 7iturcn Toesolve Bon.) 3.32, which for s = i may be 

written as 

PPire = T A Cele re Bae cist) 
where {C*} represents an arbitrary set of initial values of 
{P,}. Integration yields 

[P,} = [H.1{c#} + {0.1 3.38 
where [H,] is defined as before, {Q,} is a vector arising 
from the integration of {B,}. Since the particular solution 
is any solution which satisfies the inhomogeneous equations, 
it 18 adequate to select 

{C*} =.-0 
Hence, Eqn. 3.38 .reduce to 

Ue een b 389 
Therefore, the final solution is formed by Superimposing the 
two -solutzons, Bans sus. oeand, 3230. 


Vee meme) Gite © es 3.40 


3.1.5 Segment Stiffness Matrix 
FOR USs anon wor.40 sbecomes 
ty ;} = Poplist eer 3.41 
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"j' corresponding to each unit variable applied at ‘i' in 
the absence of any external loads. {Q,;} represents the 
variables at 'j' corresponding to zero displacements, D, and 
forces, F at 'i' in the presence of the external loads. 


Thus, Eqn. 3.41 can be expanded into 


D; Hy He D; Qu 
ate 


F, H; H,|IF;| [of 


where Qy and Qf are the displacements and forces from the 
particular solution respectively. The total matrix in Eqn. 
3.42 1s usually referred to as a ‘transfer matrix’. 
Expanding Eqn. 3.42 into two equations 
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Solving for <D, F,;> and substituting into 3.44 yields 
F; D |; 0 
ee Yo boils, “i 3.45 
vy Da-Qe OF 
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Pj Dj; Fo j 


where the coefficients of [K] represent the forces at each 
shell edge due to a unit displacement at each end white all 


other displacements are restrained. This matrix is known as. 
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the stiffness matrix. {F)} represents the fixed.end forces. 


3.1.6 Stiffness Matrix Sign Convention 

In the derivation of the element stiffness matrix and 
the fixed end stresses, the sign convention used corresponds 
to that generally used in shell theory as given in Fig. 2.2. 
As a result, the stiffness matrix will have some negative 
elements on the main diagonal. This can be corrected by 
adapting the so-called 'stiffness matrix sign convention'. 
This sign convention is shown in Fig. 3.3. It can be seen 
that the positive direction of the top normal in-place force 
N,, the top tangential shear force T,, the bottom moment M,, 
and the bottom transverse shear S, have been changed to the 


opposite direction. 


3.2 The Flexibility Approach 

The solution to the basic shell equations may be split 
into two parts, namely: the particular solution, which can 
be simplified to the membrane solution with negligible loss 
of accuracy, and the homogeneous solution which considers 
the bending stresses. This procedure is analogous to the 
flexibility method of analysis for a statically 
indeterminate structure. Program FLEXSHELL was developed 
based on this approach. To simplify the shell equations and 


limit the particular solutions, the following assumptions 


will be made: 


1, Loads are axisymmetric, i.e., 0/06 = 0, pe = 0, 
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thus, 0¢ = dd; 
The shell segment has uniform thickness; and, 
z (from Eqns. 2.5 and 2.6 is small compared with the 


Prague OL CUrVatCULe),, 1.6.4, ratze= 5, wanaers tsetse. 


Thus, the equations of equilibrium become 
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Gr eOovm te Nor ,Sing + foNe. «fotyDsn =O 35475) 
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~d(roMg) + ryMgcos¢d + Ooror = 0 3.47(c) 
don 


and the force-displacement relations become 


No = [1 (av-v ) + »_(veosé-wsing) | 3.48(a) 
Lj dd To 
No = |» (av-w) + 1 (vcose-wsing)| 3.48(b) 
r’, \d@ Lo 
Mg = D1 d (J o ) + voos¢ du] 3.48(c) 
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Me = -D] >» s(2 ot ) a cosean| 3.48(d) 
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The method of analysis is outlined as follows: 
Determine the particular solution forces and the 
deformations at the edges of the shell due to the 
applied loads; — 

Establish the flexibility matrix; 
Solve for the edge forces and moments necessary to 
restore the incompatibilities of the deformations 


between adjoining elements; 
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4, Determine the final stresses by superimposing the 
particular solution stresses and the stresses due to the 


incompatibilities. 


3.2.1 The Particular Solution 
As mentioned earlier, the particular solution is 
approximated by the membrane solution. The membrane theory 
of shells approximates the solution to Eqns. 3.47 and 3.48 
by neglecting the bending components, based on the 
assumption that the displacements due to the membrane 
stresses do not induce any appreciable bending. Thus, Eqn. 
3.47 and 3.48 reduce to two equations with two unknowns as 
shown: 
Ur oN eur Ne COSsdé i cig linbar uo 3.49(a) 
TaNeSing) ter oNe. fubol pee sn0 3.49(b) 
The in-plane forces Ng and Neg are obtained more simply 
from the vertical and normal equilibrium of the statically 
determinate shell segment under the applied loads. Since the 
Padiieof curvature £3 and, 2) vary jin ae depending on the 
typerof shell, of revolutions, so,doesethe form of, the 
membrane solution. 
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Neo= +R Peper SD eb) 
27r,Sin’7¢ 
See eCone 
N, = -R 
27SCOSAa 3,52)(a) 
Ng = *Pz2r2 3.52 ¢b) 


where R is the total vertical load, positive when directed 
toward the supports; p, is the component of the external 
load per unit area normal to the shell surface in the 
direction towards the axis of revolution. The upper and 
lower signs relate to Figs. (a) and (b) respectively, of 
Tables 3.4 and 3.6. The expression for the membrane in-plane 
forces for the spherical, cylindrical, and conical segments, 
subjected to various loading conditions shown in Tables 3.4 
boo .6 were derived from Eqns. 3.50 to G52. The solution 


due to the thermal effects were obtained from Billington(3). 


3.2.2 The Homogeneous Solution 
Consider the vertical equilibrium of a shell element, 
then 
27 EGNig SING “Helm oO gcos¢e + R = 0 
from which 
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where R is defined as before. Note that the second term is 
the membrane force which can be evaluated separately as 
shown earlier. Therefore, the homogeneous solution is 
obtained by solving the simplified shell equations (Eqns. 


3.47 to 3.48) ignoring all load terms. Thus, the homogeneous 
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solution for the meridional force is 
Ng = -Qecot¢ 3.54 
substituting this into Eqn. 3.47(b), ignoring the load term 


Pz, and using the relation, 


Yo = r2sing 3.55 
then 
No = -r2 dQg Shy oie 
r, d¢ 
Let 
U = r2Qo cow) 
V = 1_[v+aw] Sse 
Yr ag 


Eqns. 3.54 and 3.56 become 


Ng = -1_U cot? 3459 
2 

Noes iol Ou epaien’ 
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Rearranging Eqns. 3.48(a) and 3.48(b), and substituting Eqn. 


S255 yields, 
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-&_|a(No-»No) | 3.64 
d@ “Eh 
Substituting Eqns. 3.59 and 3.60 into Eqn. 3.64 yields one 


equation with U and V terms only. 


2 g20 as ‘8 (E2)* r2cot¢? at 2 a4 


hint d¢:* r,|ado\r, ren r,h d¢d|d¢ 
“Ap En COCo, 8b “= ay dhscoto) Usa" EnV 3.65 
Eek 2 h d¢ 


Substituting “Eqns. 13.57 andio3 -58)cinto Eqns) 3448c ) rand 


3248: Gd:) } 
Mo = -p[v cote + y av| 3.66 
C2 ry de 
Mo = -p| »vcot¢ + 1) av| 3.67 
Te iy d¢ 


Substituting these two equations and Eqn. 3.57 into 


Eqn. 3.47(c) yields, 
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Eqns. 3.65 and 3.68 permit a closed form solution of 
the equations of shells of revolution. The solution of these 
equations may be further simplified by applying the 
geometrical properties of each shell. 

For the cylindrical segment with r; = © and ro = rz = 


r, these equations reduce to the form (See Appendix A for 
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d*Ay + 4B*A,= 0 3.69(a) 
ds‘ 
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Biro als) 3.69(b) 
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for which the solution can be expressed in closed form as 


Au= e8*(C,cosfs + C2sinfs) + e~ 8§(C,cosBs + Cysings) 


Forethesconical segment with’ r,)= SS. Sing, r4 2.2, 90> = 
S$ tana, and @ = 2/2 - a, the closed form solution in terms 
of the Kelvin functions ber, bei, ker, kei as shown in 


detail in Appendix A is 


Q, = 1(C,ber2& + Czbeiz—& + C3ker2& + Cykeiz€é) Se Jal 
Ss 
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eae 1A ig ess Sree 
h*tan’a 
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For a spherical segment with r, = rz = a and ro = a 


Sing, Eqns. 3.65 and 3.68 become 
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Assume that the bending effects are Significant over 
only a short distance from their point of introduction. 
Thus, BQns. aoe) 4mand 3.75 Can@begreducédito one cuigrerential 


equation in terms of a single variable, for which the closed 


form solution is a function of 
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et \*{cosk¢,sindr¢} 
where A is large and dimensionless. Note that each time the 
solution is differentiated with respect to ¢, the result is 
a multiple of the large parameter i. Consequently, the 
second derivative will be two orders of A greater than the 
Solution itself and so on. Therefore, 


GeOgu hee dOgthe>s0% 
dg? d¢ 


So all lower order derivatives with respect to ¢@ may be 
neglected in the formulation of the final solution. This 
assumption was first introduced by Geckeler in 1926. Hence, 
it will be referred to as Geckeler's assumption (2,9). Thus, 


Eqns. 3.74 and 3.75 reduce to 


COs o#: EAV SEAL 
dg? 
dqd?V = -a’Qg Bist sf 
dg? 


Combining these to eliminate V, 


d*Og+ 4405 = 0 SS 
ag * 
where 
NewS al oe) are 3.79 
h? 


The final solution is 


Qo = er*(C,cosag + CzSindrg) + e- **(C3cosdkgd + Cysindr?¢) 
. Sheol) 


whnererC., Co, Cos and GC, are arbitrary constants of 
integration. The limitations of this approximation will be 


discussed in Chapter 5. 
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3.2.3 Segment Flexibility Matrix 

The>construction* of theoflexibility) matrix for the 
cylindrical, conical, and spherical segments will be 
discussed in this section. Note that for the spherical 
segment, Geckeler's assumption will be used as indicated. By 
Seminition, the flexibility matrixscoetficient, say E(, 4), 
is the deformation of the segment at i due to a unit value 
of the load applied to the segment at j. 

Only those deformations which violate continuity and 
the corresponding forces which produce these deformations 
need be identified in the formulation of the flexibility 
matrix. For axisymmetric loading, these are the horizontal 
displacement A, and the meridional rotation Ag and the 
corresponding forces H and Mg at each discontinuous edge of 
the shell. Combining Eqn. 3.76 and 3.67 and again for the 
sphere neglecting the lower order differentials with respect 
to ¢, the expression for the meridional moment becomes 


Mg = =D_ d*Qg 3.81 
Eha d¢? 


and from the geometry of the shell, the horizontal force, 
Ha, ican be erpneedes as a function of the meridional force 
Ng. Thus, for a spherical segment 
H = Ngcos¢ 3.02 
Consistent with the sign conventions shown in Fig. 3.4, 
expressions for the moment and horizontal force at each 
shell edge may be expressed in terms of the homogeneous 


solution, as shown in matrix form below for the spherical 


segment (Eqn. 3.74) 
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Let 

@, = efcosk¢ 6, = e**(cosr¢ + sindd) 

2 = esindd 62 = e*(cosrk¢ - sindrd) 

3 = e Scosr¢ 6, = e-A*4(cosrd + sindd) 

dy = e A*Ssinde 6, = e -A>*(cosd¢@ - sindd) 

then 

Hy =\/sinao. 0 0 0 1 92 63 H/C, 
M6 Uo aye 0 0 ei Mel Oe CR rele 
Hi] 0 0 1/sine 0 || 61 64 04 oi| |e 
Md 0 0 Ge a/2n1 1-61 64-6) 61) 1c, 


Or simply, 

EV ir t= ThT PRT Be} 
Multiplying matrices [T,] and [Tz] simplifies to 

(Viiee= OORT GC} a8 
Similarly, expressions for the deformations at each shell 
edge may be expressed in terms of the homogeneous solution 


as follows: 
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Ag 0 0 O 2d?! |-94 1 o4 -O4]]Cu 


Or simply, 
PA ee=e tr. iT, 116s 
Multiplying matrices [T;] and [T,] yields 
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Combining Eqn. 3.83 and 3.84 yields 
iA} 
{A} 


where [F] is the segment flexibility matrix, such that 


BRA rr | -asty } 3.85(a) 


[F]{v} 3.85(b) 


fel) Ss [Tai lra)-+ 3.85(c) 
Similarly, the flexibility matrices for the cylindrical and 
conical segments are constructed using the homogeneous ~ 
solutions, Eqns. 3.70 and 3.71 respectively, as shown in 
detail in Appendix B. 

The base segment is considered to be a circular plate 
Supported on a Winkler type foundation, whose stiffness is 
expressed as the subgrade modulus, k (4). The segment 
flexibility matrix was developed in the same manner as the 
spherical, cylindrical, and conical segments, based on the 


asymptotic solution to the fourth order plate equation. 
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TABLE 3.1 Coefficients of Matrices B,, Bg and Load Vector B3 in Eqn. 
3.19 


Matrix Bo 


Vector B3 
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TABLE 3.2 Coefficients of Matrix A 
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TABLE 3.3 (a) Coeffiecients of Matrix By in-Equ~3627 
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Table 3.4 MEMBRANE SOLUTION FOR A SPHERICAL 
SEGMENT 


MERIDIAN 
— UPPER AND LOWER SIGNS 
PARALLEL Z-AXIS RELATE TO ae (a) AND (b) 
CIRCLE e RESPECTIVEL 
N-«— — Wi & 
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Table 3.4 (cont'd) 


LOAD CASE IN-PLANE FORCES & DEFORMATIONS 
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Table 3.5 MEMBRANE SOLUTION FOR A CYLINDRICAL 
SEGMENT | 


LOAD CASE 
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Table 3.5 (cont'd) 


LOAD CASE IN-PLANE FORCES & DEFORMATION 
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Table 3.6 MEMBRANE SOLUTION FOR A CONICAL 
SEGMENT 


—- UPPER AND LOWER SIGNS 
RELATE TO Figs (a) AND (b) 
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Table 3.6 (cont'd) 


LOAD CASE IN-PLANE FORCES & DEFORMATIONS 
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Figure 3.1 EFFECTIVE SHEARING FORCES EXPRESSED AS A 
FUNCTION OF THE IN-PLANE SHEAR AND THE 
TWISTING MOMENT 


TANGENT TO THE POINT 
BEFORE DEFORMATION 


TANGENT TO THE POINT a 
AFTER DEFORMATION Jr 


Figure 3.2(a) MERIDIONAL ROTATION B DUE 
TO DISPLACEMENT 7 


Figure 3.2(b) MERIDIONAL ROTATION B DUE 
TO DISPLACEMENT a 
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Figure 3.3 SASHELL STIFFNESS MATRIX 
SIGN CONVENTION 


54 


SHELL EDGE =i 


- SHELL EDGE j 
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Figure 3.4¢ FLEXSHELL FLEXIBILITY MATRIX 
SIGN CONVENTION 
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4, FLEXSHELL FORMULATION 

Program FLEXSHELL analyzes a multi-shell structure 
based on the flexibility approach described in Section 3.2. 
The program listed in Appendix C is a modification of an 
earlier version developed by Murray, et al(4) for the 
analysis of the Gentilly type containment structure. The 
capability of the program has since been increased for a 
wider variety of multi-shell problems. 

The 'long' sphere was replaced by a 'short' spherical 
segment. And three additional segments were introduced such 
as: the 'short' conical segment, the 'short' inverted 
Spherical and conical segments.(Fig. 4.1) Furthermore, two 
load cases were added, namely: the liquid pressure loading 
and the snow load, which is a uniform pressure over a 
horizontal projection of the shell segment. Consequently, a 
Significant portion of the original version of the program 
was recoded. Specifically, 

1. The membrane solution forces developed in Chapter 3, 
shown in detail in Tables 3.4 to 3.6, were coded 
directly into the program. The upper and lower~signs 
relate to the dome and inverted dome configurations. 

2. The construction of the flexibility matrix was coded 
from the product 

Ere) Saray 
whose coefficients were coded directly into the program. 
Provisions were made for the dome and the inverted dome 


configurations and also for the closed spherical and 
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conical segments. 

The calculation of the particular solution displacements 
were coded uSing the expressions for A, and 4g, shown in 
Tables 3.4 to 3.6, using the membrane solution forces 
already coded into the program (Step 1). 

The special homogeneous solutions (4), which considers 
the effect of the vertical edge load was replaced. 
Finally, the calculation of the final stress resultants 
and displacements were coded to conform to the 'short' 
segments that were introduced. 

The logic flow of program FLEXSHELL is as follows: 
Define the segment connectivities; 

Satisfy the rigid body motion requirement by evaluating 
the effect of the vertical load components on the 
segment; 

Evaluate the joint eccentricity effects; 
Establish the segment flexibility matrix; 

Solve the compatibility equations; 

Solve for the final stress resultant values by 
superimposing the particular solution stresses with the 


stresses due to bending. 


4.1 Definitions and Notations 


Segments are defined with reference to a coordinate 


starting at the line of symmetry at the apex of a structure, 


andutraversingsthe midsurface or ther shell segment in ja 


counter clockwise sense, until again reaching the symmetry 
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line at the base of the structure. The coordinate for 

branches which do not fall in this primary circuit may be 

defined in the same manner, starting at the free edge, 
increasing in a counter clockwise sense. Therefore, with the 
exception of the last segment in the primary circuit, the 

"bottom' of each segment is always supported by the 'top' of 

the adjacent segment. For reasons which will be explained 

later, the segments must be numbered sequentially in such 
manner that any segment always has a higher number than any 
of the segments which it supports. 

Consider a shell segment cut by a vertical plane shown 
in Fig. 4.1, the sign conventions consistent throughout the 
program are as follows: 

1. Moments and rotations are positive as shown in the 
figure. 

2. Horizontal forces, displacements, and eccentricities are 
positive in the direction towards the line of symmetry, 
also known as the axis of revolution. 

3. Vertical forces are positive downward. 

4, For base segments, vertical displacements are positive 
downward, whereas, vertical eccentricities are positive 
upwards. 

Note that all forces and moments are expressed per unit 


length. 
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4,2 Connectivity Matrix 

The connectivity matrix is established by satisfying 
the geometric compatibility requirements between adjacent 
segments. This is accomplished by forming the algebraic 
summation of the horizontal displacements and meridional 
rotations at adjacent edges of the shell segments. To 
express these equations in matrix form, it is necessary to 
number the segments as described earlier, to ensure the 
consistency in the order of assembly of the end 
deformations. Furthermore, associated with each segment is a 
flag indicating the presence or absence of a connection at 
the 'top' and 'bottom' of the segment, input as IR and UR, 
respectively. The connection between segments is specified 
by POCO eL.-l) and TUCOUl"; 2) Whien 1S the number Of elie “bop ) 
segment and the adjacent ‘'bottom' segment, respectively. The 
compatibility equations expressed in matrix form is 

[A] {A}, = {0} 4.1 

where [A] is the Boolean connectivity matrix expressing the 
compatibility requirements between segment deformations (4); 


fA}, is the total segment deformation vector. 


4.3 Vertical Edge Load 

This section demonstrates how the rigid body motion of 
the shell structure which have been ignored up to this point 
is taken into account. Loads from the 'top" segment may be 
transmitted to the segment below it as a vertical edge load 
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which can carry this load by membrane action alone, for the 
case of the spherical and conical shells, a horizontal force 
H, must be added vectorially, so that a resultant force Ng 
is formed (1,6,14). This horizontal force must be 
compensated later by subtracting this value from the real 


horizontal loads PSF(N,1) and PSF(N,3) acting on segment N. 


4.4 Shell Eccentricity 

Since segments at a joint may not always end at the 
Same point, a horizontal segment eccentricity may be 
specified in the input data. This results in the 
eccentricity of the edge horizontal and vertical loads, 
which in turn produces a moment which must be added to the 
existing moments PSF(N,2) and PSF(N,4) at the edges of 
segment N. This moment is automatically calculated in the 


program. 


4.5 The Particular Solution 

The particular solution is approximated by the membrane 
solution. The computation of the membrane in-plane forces Ng 
and No, for the spherical and conical segments are 
incorporated into. function: subprograms: FN1/* FN2,. FN3),. and 
FN4, respectively. The equations used in these subprograms 
aretfounduan’ Tables se40and S2l6é-ewThersotutioni torathe 
cylindrical segment, found in Table 3.5, is simple enough, 
that a separate subroutine is not necessary. The particular 


solution displacements PSD are obtained from evaluating the 
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equations for A, and Ag found at the end of Tables 3.4 to 
3.6. These computations are incorporated into subroutines 
PCYLIN, PDOME, and PCONE, respectively. The particular 
solutions for the base segment derived in (4) are 


incorporated into subroutine PBASE. 


4.6 The Flexibility Matrix 

As derived earlier, the flexibility matrix for a shell 
of revolution may be expressed as follows 

Pear ACh aries 

These matrix operation is performed by subroutines CYLIN, 
DOME, CONE, and BASE, for the cylindrical, spherical, 
conical, and base segments respectively. 

The first step is to initialize the coefficients of 
[TA] and [TT]. For subroutine CONE, this necessitates the 
use of another subroutine MMKEL2, which computes the Kelvin 
functions of order 2 using published recurrence formulas 
(10). Subroutine MMKEL2 in turn calls up a system-dependent 
Subroutine which evaluates the Kelvin functions of order 
zero and one and their derivatives. Secondly, a check is 
made if the segment is inverted or not. By definition, an 
inverted spherical or conical segment is that which forms a 
cup-like shape as shown in Figs. (b) of Tables 3.4 and 3.6. 
If the segment is inverted, subroutine ROWEX is called. This 
performs row interchanges in the [TA] and [TT] to conform to 
the inverted configuration. Furthermore, a check is made 


whether the segment is a closed spherical or conical dome. 
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If so, the four by four flexibility matrix degenerates into 
a two by two matrix. The next step is to invert the [TT] 
matrix which is performed by subroutine TTINV which is 
capable of inverting a four by four or a degenerated two by 
EWO Matrix. Finally, the matrix multiplication 

EATER Tse 


is performed, thus forming the flexibility matrix. 


4.7 Matrix Formulation of the Solution Procedure 

Let [F],; be the flexibility matrix of segment i, then 

Peomesan. 3.85(b)’, 

tA ae ve 4.2 
Similarly, for the entire structure, the equations are 

Ay = IV 4,3 
where the end displacements {A},, end forces {V},, and 
flexibility matrix [F], of element i are assembled into the 
global matrices {A}, {Vv}, and [F], respectively, in the 
order consistent with the sequence of segment numbering. 

The particular solution displacements and the vertical 
edge load displacements {6} in the corresponding order as 
{A}. The total displacement vector is 

CASH eae tO 4.4 
Substituting Bon. 4. 4winto 451, 

PAI CUA Seay) = 103 4.5 
and multiplying Eqn. 4.3 by [A], 

[AJ(F]iv} = [A] {4} 4.6 


Let {q} be a set of relative displacements in terms of the 
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homogeneous solution {A} such that 

{iq} = [A] {4} 4.7 

From a general theorem in structural analysis (13), if 

a set of forces {V} is associated with a set of 
displacements {v}, and if in another coordinate system, the 
Same set of forces may be described as {U}, and their 
associated displacements as {u}, then the work done in the 
two systems must be identical when undergoing equivalent 
displacements, i.e., 

<u Uy) tv 4.8 
Similarly, 

= 0s = se ive, 4.9 
where {V} are the forces associated with displacements {A} 
and {Q} are the redundant forces associated with the 
relative displacements {q}. Substituting the transpose of 
Bon. 4.7 into’ 4,9) 

<ASER | fO} =) <A> iG 4.10(a) 

<A>({A]'{Q} - {v}) = 0 4.10(b) 


since Eqn. 4.9 must be true for all <A>, Eqn. 4.10(b) 


becomes 
{v} = [a]' {a}, Nel 
Substituting eqn. 4.11 into 4.6"yieids 
[A][FICA]' {Q} = -[a] {6} 4.12(a) 
[F1{Q} = {qo} 4.12(b) 


where the structure flexibility matrix and structure 
particular solution displacements, respectively, are 


Pee Gan) LES Pade 4.13 
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igo} = -[A] {8} 4,14 
The set of simultaneous equations (Eqns. 4.12(b)), can 
then be solved for the redundants Q, by Gauss elimination. 
Once evaluated, the value of the redundants may be 
back-substituted into Eqn. 4.11, to find the edge forces V, 
whieh @Raturnecan be -substituted into Banw.4.3, to find: the 
displacements A. The final solution can then be obtained by 


Superimposing these values on the particular solution. 
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Figure 4.1 SUBSCRIPTING OF SEGMENT 
ARRAYS FOR FLEXSHELL 
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5. EVALUATION OF THE FLEXIBILITY APPROACH 

The flexibility and stiffness methods are two basic 
approaches in the analysis of statically indeterminate 
Structures. For a given problem, both methods will give 
identical solutions. Any differences which may be observed 
are due to the approximations used in the formulations, and 
are not inherent in the methods themselves. Basic 
limitations such as axisymmetry of the loads and geometry, 
if introduced, will limit the scope of the applications of 
each method. 

The Simplifications used in formulating the equations 
on which program SASHELL is based are consistent with the 
elastic theory. Hence, any errors in the solution are due to 
the manner in which the equations are solved. In SASHELL, 
these may result from the loss of stability of the numerical 
integration technique used, and, for non-axisymmetric 
loading, the rate of convergence of the Fourier expansions 
could be a factor. It should also be noted that the elastic 
theory is based on an element that has two edge boundaries 
along®the meridian, So insthesformulationwof the stiffness 
matrix, four boundary conditions must be imposed at each 
edge of the shell in order to evaluate the eight constants 
efsintegrabion mThus i) theactosedpespher 1calieor ec emi call 
segment which has only one edge results in a Singular 
stiffness matrix. This problem may be overcome with the 


introduction of a small hole at the apex, say ao = 0.6°. 
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On the other hand, the solutions obtained from program 
FLEXSHELL are based on closed form solutions. The degree of 
Simplifications introduced in obtaining these solutions 
depend on the shell type. With cylindrical shells, no 
Simplifications are necessary. Hence, the solution should be 
identical to that from program SASHELL. With conical shells, 
some error may be introduced depending on the semi-vertex 
angle, a due to difficulties in evaluating the Kelvin 
functions and their derivatives. Using Geckeler's 
approximation in obtaining the closed form solution for the 
spherical shell will likely result in greater error, 
particularly as the values of A or a become small. 

Closed form solutions for each segment were obtained 
for 'short' segments, that is, two constants of integration 
are evaluated at each edge of the shell. However, unlike 
program SASHELL, the analysis of closed, spherical or 
conical segments offers no difficulties. In such case, only 
the constants of integration corresponding to the end 
furthest from the apex are evaluated. Thus, imposing 
boundary conditions at the apex are not required. 

In order to evaluate the reliability of the two methods 
used in the programs, both programs were run on identical 
problems and the results are compared. Comments as to which 


is considered to be more reliable for a specific problem are 


presented. 
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5.1 Cylindrical Segment 

Fig. 5.1 illustrates the distribution of the bending 
moment and circumferential force along the length of the 
cylinder, under a hydrostatic pressure with y = 62.4 pcf, 
obtained from programs FLEXSHELL and SASHELL. The meridional 
force N, is equal to zero in this case. As expected, the two 
solutions are identical. Similarly, upon investigating the 
solutions for the other load cases, dead load, snow load, 
and uniform pressure or prestressing, both programs yield 


id@entical results? 


5.2 Conical Segment 

There are phySical limits that must be imposed on the 
semi-vertex angle a for conical segments. AS a approaches 
zero, the cone degenerates into a line and no shell action 
is possible. On the other hand, as a approaches 90°, the 
cone becomes a circular plate. 

With the formulation of SASHELL, the latter case 
presents no problem since the basic shell equations reduce 
to the plate equation when a = 90°. However, in FLEXSHELL, 
ibeisonecessary*towimpose a limit onethesrange ofuvaluestof 


£ for which the Kelvin functions and their derivatives are 


evaluated. 
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approaches ©. Thus, this dimensionless parameter used in the 
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evaluation of the Kelvin functions and their derivatives are 
limited as follows: 

Ovetéiose 149.40 
Gonsequently, a limit on the vrange*of-valuesvof ais 
imposed. For instance, when s/h = 500, a must be greater 
than 26°, and when s/h = 100, a must be greater than 5.5°. 
Note that a can be very close to, but not equal to 90°, say 
SIS” : 

Fig. 5.2 illustrates the distribution of the in-plane 
forces Ng and Ne and the meridional bending moment Mg along 
the conical segment, under snow load, p = 1 ksf, according 
to programs FLEXSHELL and SASHELL. As anticipated, both 
programs yield identical results since no simplifications 
were made in the formulation of the closed form solution for 
this segment. The solutions for the other load cases are 


also identical. 


5.3 Spherical Segment 

Fig. 5.3 ibbustratesetherdtstmnbuttoneot ithe munaplane 
forces Ng and Ne and the meridional bending moment Mg along 
the spherical segment for « = 10° and a = 80°, under dead 
Koadtwibh ys= elbOdpcipyaccording Lt onpnograns INLEXSHELEwand 
SASHELL. It is observed that there is a greater discrepancy 
between the solutions for small values of the angle a, say 
0%, sthanifom bargetvallues ioflapasay abOvonThas Gbservation 
is confirmed with the investigation of the solutions for 


several values of a with ao = 0°, aS Shown in Fig. 5.4. The 
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solutions for different values of a/h were also compared, 
and a greater discrepancy is observed for small a/h values. 
These observations may be explained as follows. As the angle 
a become small, the bending effects become more significant 
over a large portion of the segment, Clearly, this violates 
the basis of Geckeler's assumption; hence, the approximation 
for the spherical segment becomes less accurate. For the 
Same reason, a greater discrepancy is observed for small 
watuestot ({asag)with ag #50°)"as, ivbustratedtinsrigre 525. 
Note that no discrepancies are observed for a = 90°. Since 
Geckeler's assumption requires that the dimensionless 
parameter A, which appears in the closed form solution for 
the spherical segment, must be large, that is, a/h is large, 
the approximation becomes less accurate for small values of 
a/h. In order to be able to predict the discrepancy between 
the two solutions for a specific set of geometry and 
material property, Figs. 5.4 and 5.5 were combined to 
develop Figs. 5.6 to 5.8 for various load cases. 

Fig. 5.6 compares the meridional forces in a spherical 
segment under various load cases as given by the two 
Gomputer programs. AS illustrated in the i1gqure; ihe 
solutions obtained from FLEXSHELL show excellent agreement 
with those obtained from SASHELL, the maximum difference 
being 1/200 of 1 percent. For the liquid pressure loading, 
the solutions become identical as a approaches 90°, and 
smaller discrepancies ar observed with increasing values of 
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Solutions differ by less than 5 percent. 

Fig. 5.7 compares the circumferential forces on a 
Spherical segment under various loads. Similarly, the 
solutions become identical as a approaches 90°, and smaller 
discrepancies are observed for high values of A(a-ag). The 
discrepancies between the solutions for A(a-a,) = 10 due to 
dead load and snow load and due to uniform pressure are less 
than 5 and 10 percent respectively. For the liquid pressure, 
a greater discrepancy is observed, up to 10 and 5 percent 
discrepancy is observed for A(e-a,) = 43 and A(a-ao) = 60 
respectively. 

Fig. 5.8 compares the meridional bending moment in the 
Spherical segment under various loads. As with the in-plane 
forces, the bending moments from both programs become 
identical as a approaches 90°, and smaller discrepancies are 
observed with increasing values of A(a-ao). A discrepancy of 
less than 5 percent for the uniform pressure, dead load, and 
snow load is observed. However, a discrepancy of less than 
10 percent is observed for any @ with A(a-a.) = 25, and 5 


percent for A(a-ao) = 45. 


5.4 Application of Program FLEXSHELL 

Having investigated program FLEXSHELL for simple shell 
segments, a multi-shell structure consisting of a 
combination of cylindrical, spherical, and conical segments 
was analyzed to demonstrate the capabilities of the program. 


Again, the results obtained from FLEXSHELL are compared with 
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those from SASHELL. The same example problem is used in the 
FLEXSHELL user's manual (Appendix C) to illustrate the input 
anc outpuc tules..Fig. 5.9 illustrates, anmantze tank 
consisting of cylindrical, spherical, conical and base 
segments, and ring beams which are modelled as cylindrical 
segments. Intze tanks are mainly used for water storage, and 
they are typically constructed as prestressed concrete. The 
material properties which were used are: 


E OSS 0S ex 02 ast 


poz Ueror 

y #150 mpc & 
The base of the tank is considered to be fully fixed so the 
base segment was given a high modulus of elasticity, 

E(Dase) en | sO er l0r cs pet 

The results of the analysis of the Intze tank under 

dead load, given by both programs are shown graphically in 
Figs. 5410 to 5.72. It is vapparent. from” these™iugures that 
the solutions show excellent agreement between the in-plane 
forces Ng and No, and meridional pendie momentMg, for each 
segment. Some discrepancy is observed in the solution for 
both spherical seqments in the vicinity of the apex. This 
may be introduced by the singularity formed at this location 
in program SASHELL, or due to inaccuracies introduced by 
uSing Geckeler's assumption for small values of the angle a 
in FLEXSHELL. Nevertheless, both solutions show good 


correlation and the same is anticipated for the other load 


cases. 
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Figure 5.1 CIRCUMFERENTIAL FORCE & MOMENT 
ALONG A CYLINRICAL SEGMENT UNDER 
HYDROSTATIC LOAD, ¥= 62, 4pcf. 
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Figure 5.2 IN- PLANE FORCES & MOMENTS 
ALONG A CONICAL SEGMENT UNDER 
SNOWLOAD, p =! ksf. 
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Figure 5.3 IN- PLANE FORCES & MOMENTS 
ALONG A SPHERICAL SEGMENT 
UNDER A DEADLOAD, 
¥ =150 pcf 
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Figure 5.4 COMPARISON OF THE STRESS RESULTANTS FOR 
THE SPHERICAL SEGMENT IN TERMS OF a/h 
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Figure 5.5 COMPARISON OF THE STRESS RESULTANTS FOR THE 
SPHERICAL SEGMENT IN TERMS OF (of - at,) 3 
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Figure 5.6 COMPARISON OF THE MERIDIONAL FORCE, No 
FOR THE SPHERICAL SEGMENT UNDER 
VARIOUS LOADS 
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Figure 5.7 COMPARISON OF THE CIRCUMFERENTIAL FORCE, 
Ng FOR THE SPHERICAL SEGMENT UNDER 


VARIOUS LOADS 
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Figure 5.8 COMPARISON OF THE BENDING MOMENT, M, FOR 
THE SPHERICAL SEGMENT UNDER: 
VARIOUS LOADS 
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Figure 5.10 
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6. SUMMARY AND CONCLUSIONS 

The two basic approaches in the elastic analysis of 
multi-shell structures were discussed. With the stiffness 
method, the stiffness matrix relating the forces and 
deformations at the edge of each shell segment are formed by 
numerically integrating the basic shell equations which were 
expanded into a Fourier series. These element stiffness 
matrices are then superposed to form the structural 
stiffness matrix from which the segment edge deformations 
are computed. With the flexibility method, the flexibility 
matrix is formed from the closed form solutions derived in 
Chapter 3. The particular solution is approximated by the 
membrane solution also derived in Chapter 3 and shown on 
Tables 3.4 to 3.6. A computer program, FLEXSHELL was 
developed based on the flexibility approach. The results are 
then compared with the results from program SASHELL 
developed by Shazly (5) based on the stiffness approach. 

Individual segments under various load cases were 
investigated. It was found that both programs yield 
1aentical solutions for the cylindrical, and conical 
segments. Since Geckeler's assumption was used in the 
formulation of the solution for the spherical segment in 
FLEXSHELL, a discrepancy between the solutions was 
anticipated and observed to be a function of A(a-ao) and the 


angle a, where X is the dimensionless parameter in terms of 


ay. 
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The Intze tank problem discussed in Chapter 5 was 
selected to demonstrate the capabilities of program 
FLEXSHELL to analyze an axisymmetric segmented shell 
structure. Overall, program FLEXSHELL showed excellent 
agreement with SASHELL. The main advantage to using 
FLEXSHELL is that input is simple. 

Therefore, it may be concluded that FLEXSHELL is a 
simple, effective tool for the analysis of a wide variety of 
axisymmetric multi-shell structures. 

Further development is possible, with the addition of 


more load cases and more shell types. 
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APPENDIX A-DETAILED DERIVATION OF THE HOMOGENEOUS SOLUTION 


For a cylindrical segment with geometrical properties 


as follows, 


Wo, Sl 9. = tL Ng = N, as = r,;q_ 
Ges a pd = Ds d¢ ds 
G2 = 1/2 


Eqns. 3.64 and 3.67 derived earlier become 


reaeU = BBtv Al 
ds? 

G-d2V = Tee A.2 
ds? D 


Combining these equations and using eqn. 3.57 to eliminate 


U, 
d*Ay + 4B*Ay= 0 Aw3 
ds* 
where 
B* = 3(1-vp2) 7 Neat 
sae as 


forswhich the final solution’ is 
Ay= e8*(C,cosBs + C2sinBs) + e° 8*(C3cosBs + Cysin§gs) 
A.5 


The geometrical properties of a conical segment is 


ro = S§ Sina @ =t17 2) -* @ ae sd 
1 Ne Ng = N, ) ds 
r2 = § tana Pe) =P; 


Substituting these relations into Eqns. 3.64 and 3.67 yield 


eoasUatedU tana = U tana (= BhV A.6 
as; ds l 2 

Oe Vee dy stants y Cdl Om ea ee) 
as ds C2 D 


and Eqns. 3.56 become 
U = s Q,tana A.8 
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Sede (GO + SOs) re sSO ea ennVvect a Ae. 
ds? ds Ss 
SAE VT MAE Vener ESO), ) 10 
ds? ds Ss D 
These equations can be solved to form a fourth order 
equation in terms of a single variable(7). However, an 
alternative approach 1S possible by introducing a linear 
differential operator as follows(1): 
Pt l= tsiid? edie id Jocks O! 11 
ds? ds S} 
thus, eqns. A.9 and A.10 become, 
ESO.) =s#8hV cota 12 
L(v) = -(sQ,) 13 
D 
Operating on Eqn. A.12, and substituting back into Eqn. 13 
yields, 
LL(sQ,) + A*(sQ.) = 0 14 
where 
ee es a er) 
Meebelie & 
This may be written in either of the following forms, 
BEbCsOe y+ pelesepn ls “ih LO SOei EN OSes ibs) 
BEoCeO.) = LN USO, at TALE So el CS Onan le sae pate 
which show that the solutions of the two second-order 
equations are 
CSO) Gein SO. es =e Sal 
Expanding this equation yields, 
SudclsO.)) } aso sO. yates. a0 A.18(a,b) 
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The solution to Eqns. A.18 is complex, and it will be enough 
to solve one of the equations, and then use the real and 
imaginary parts of this solution separately as the solution 
of a fourth-order equation mentioned earlier. Introducing a 
new variable, 

N= 2AV 15, A.19 
Eqn. A.18(a) become, | 


d?(s0.) +2? a€(s9,)\ + (' - 4 ) (sQ.) =e 0 A.20 
dn? n dn Be 


The solution of this equation consists of Bessel functions 
of the second kind. 
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Let & = 2A/Ys, then rewriting Eqns. A.21 in terms of the 


Kelvin functions of order zero yield 


J2(n) = 2bei'f-beré + Vpeas as) A.23(a) 
E E 
HS‘)? (n) = 2 2ker'£+keié 42 BINED) A.23(b) 
n & uy 


These two functions are independent solutions of Eqn. A.18, 
and their real and imaginary parts separately will satisfy 

the fourth-order equation formed by combining Eqns. A.9 and 
A.10. The general solution for a conical shell is 


QO, = ia, (beré-2bei'£) + (ONS aa 1) 
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Using the recurrence formulas(10) for the Kelvin 


functions Eqn. A.24 can be rewritten as follows: 
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APPENDIX B-CONSTRUCTION OF THE SEGMENT FLEXIBILITY MATRIX 


In general, the flexibility matrix of a shell segment 
is of the form 
CE =) fTad [red = Bel 
The [TA] and [TT] matrices are a function of the geometrical 
and Material properties of the shell segment. 
Based on the geometrical properties, the [TA] and [TT] 
matrices for the cylindrical segment, derived in a similar 


manner as for the spherical segment in Chapter 3, is as 


follows 
bee 
pees (1-v?) 
bene 
D= Eh? 
eiGten 

and 
¢, = e8cosfs 6, = e®* (cosBs + sings) 
@, = e8sings 6, = e®*(cosfs - sings) 
$3, = e- 8cosfs 6, = e 8*(cosfs + sings) 
¢, = e- 8*sings 6, = e- 8§(cosfs - sins) 
H! 2DB? 0 0 0 }i=4 1 1 bo Meee 
Mo 0 FABY 03 2 0 0 0 am 0 1 C2 
ik ghuegs = domereney 0) |) eles eh. oe anno. 
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The [TA] and [TT] matrices for the conical segment is 


Ss follows -(8); 
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Inverted shell 


PrAshies B.5 


An inverted cone or sphere is shown on Fig. (b) of 
Tables 3.5 and 3.6 respectively. Note that in the above 
derivations, i and j, relates to. the 'top* and ‘bottom’ of a 
shell segment. So, for an inverted shell, the top becomes 
the bottom and vice versa. Thus, to find the flexibility 
matrix, it is a simple matter of interchanging rows one with 


three, .and rows two with four, of, matrices [TA] seq. [TT]. 
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APPENDIX C-FLEXSHELL USER'S MANUAL 


Using the flexibility method of analysis, program 
FLEXSHELL computes the in-plane forces, bending moments, and 
horizontal displacements for an axisymmetrically loaded 
Shell structures due to various loads. 

The program is capable of analyzing six types of shells 
of revolution of uniform thickness. These are cylinders, 
spheres, inverted spheres, cones, inverted cones, and base 
slabs on an elastic foundation. Seven axisymmetric loading 
cases are available. These are self weight, uniform 
pressure, prestressing, snow load (a uniform vertical load 
over a horizontal projection), hydrostatic load, uniform 
temperature change, and temperature gradient through the 
shell thickness. 

The input to FLEXSHELL consists of multiple lines which 
may be lines in a datafile or a set of punched data cards. 
There are six input card types. Certain card types may be 
repeated as necessary. 

A typical explanation of a card type consists of the 
card type number, a descriptive name indicating the nature 
of the data, the format used, and the number of cards of 
that type) required nthisoised obfowedabytthe yarmable names, 
in bold type, followed by the definitions of these variable 
names, and the options available ,if any, for the input 
variables. Throughout the input all units have to be 
consistent. The input and output files for the Intze tank 


problem discussed. in Chapter 5 and the program listing are 
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given in the latter part of the Appendix. 


ee 


TITLE card Format 10A8 
Any identifier string up to 80 characters. 

CONTROL card Format 214 
NSEG IPRINT 


NSEG = number of shell segments in the structure. 


IPRINT = print control character 


0 - echos input data and prints final results only; 

loz lppints Bil routputeneludings the connectivity 
matrix, PSF and PBF arrays, and the element and 
Structure flexibility matrices. (used for 
checking purposes only) 


SEGMENT DATA card Format 514,2F10.4 


One card per segment required. Note that the segments 
must be numbered sequentially in such manner that any 
segment always has a higher number than any of the 
segments which it supports. 


Port sir (1) sDRC2)INDLV, BEC ine GGle 2) 
I = segment number 
IT = segment type 


1. cy langes 

2.7 sphere 

3 - base on elastic foundation 
4 - cone 

5 - inverted sphere 

6 - inverted cone 


IR(1) = top connectivity flag for segment 1 


OT -tOp use not cOonnectea CoO vanotner segment, 
1°- top 1s connected to another segment. 


IR(2) = bottom connectivity flag for segment I 


O° ="bottom=1s not connected to another’ segment; 

1 - bottom ys connected to another segment; 

-1 - bottom is connected to another segment with a 
pure hinge. 


THs 
" a 7 
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NDIV = number of divisions for segment I at which stress 
resultants are to be computed and printed. (max = 100) 


EC(I,1) = eccentricity of joint connection at the top of 
the segment in feet. 


EC(I,2) = eccentricity of joint connection at the bottom 
of the segment in feet. 


NOTE: 

When two shell segments are connected at a given 
elevation but have different midsurface radii, a 
horizontal eccentricity equal to the differences in 
horizontal radii to the midsurfaces will result. This 
eccentricity can be applied to either shell segment and 
1s positive when directed inwards. For the eccentricity 
between a spherical and cylindrical segment, EITHER of 
the following entries is permissible. 


EC 61,22 )\4= 40 EGG, 2) e> =oad 
EC (271). =ecno BCC2, (rea 0 
CONNECTIVITY specification card Format 214 


Specifies the connection between segments. Requires 
(NSEG - 1) cards. 


IDCO(1) IDCO(2) 
IDCO(1) = number of top segment. 


IDCO(2) = number of segment to which the top segment is 
connected. 


NOTE: 
Where three shell segments intersect at the same 


elevation, two connectivity specification cards are 
required. For the Intze tank shown in Fig. 5.9% the 
entimesfaressear, on the first card, andtGre4qron ithe 
following card. Each segment number appears precisely 
once in IDCO(1) and these numbers must be arranged 
consecutively in increasing order, starting with segment 
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1 and ending with segment NSEG. 


SEGMENT PROPERTIES card 


I T R H HO E PR ALPHA UW 
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Per Oye Un 


I = segment number 

T = segment thickness, feet 

R = radius of the parallel circles for cylinders and 
spheres, feet; or, 
= subgrade coefficient for base segments; or, 
= semi-vertex angle of cone in degrees. 

Hoe lengen sin feet. for cyl wnder-s-0or, 


= total angle in degrees from the axis of revolution 


to the outer edge for a sphere; or, 


= outer radius of a circular base ring slab in feet; 


OL, 


= distance from the apex of cone to the large end in 


feet. 


HO) =-.0/50 or, blank for cylinder-‘or, 


= angle in degrees from the axis of revolution to 


the inner edge for a sphere; or, 


= inn@rs radius Ofwas CaAnculam basemings skab: or, 
= distance from the apex of cone to inner edge in 


feet; or, 
E = Young's modulus for segment, psf 


PR = Poisson's ratio for segment 


ALPHA = coefficient of thermal expansion 


UW = unit weight of material for segment, pcf 


LOAD TYPE card 


One card per segment is required. 


Format 214,7F10.0 


- 
wi caete ae sla siecle 
us / _— , 7 

Mea | Anes dnemgea or 

is ae ; 7h : 

oan: Na dawiaabis: a nega » bie 
. So i 5 a : o 
ia. +74 OO mniarep) Eg Peete e e <3. 7 ‘aeet x RA ~. 

ag te Cian 387 wtiQe 


20 :8§eehpet gehd 367) ‘Sota io ttieeat, Aheigaus =. 
jaee Geb . Tr) 4105 Fo, ahere. a th well ia ig = a 
wh Bey ; / 
ue ity ? ez ven bith Apphes as * 
2teeh yb edge Laged -« 
1 ye See e388 eobe setio arts of \ 
i 6845 SHOT Pies ca = GORE au ibaz. yest. - 
ne ; ¥ : ¢ - ah a. .. an - tal 
G1 ° Shs, OF GAGS te ‘on ade waad—wan ates 
iy ae nis. rad aw sae) 

wae ye 

| ons) 420 coi steve 
AC i SUPOVEaR tore 2 dk Bs 24899 ee } 


5 SIs Teli 
> ae et Te Ye s 
= ie beer ta 
ea, ot 


ie 


t 
‘om 


101 


I IP PV WHT PSF(1) PSF(2) PSF(3) PSF(4) PSF(5) PSF(6) 


I = segment number 


IP = load type parameter 


uniform pressure 

self weight 

prestress loading 

uniform temperature change across section 
temperature gradient across section 
uniformly distributed load over a horizontal 
projection, or snow load 

liquid pressure 


NOTE: 


A hydrostatic load applied to the base segment 


is Simulated by using a uniform pressure equal to 
the product of the liquid weight density and the 
height of the water above the base. 


PV = value of the applied load, depending on the type of 
load. 


If 


BE 


Lf 


if 


LE 


Pe 


IP=1, PV is the magnitude of uniform pressure in 
psf. 

Positive for internally directed pressure and 
negative for externally directed pressure. For a 
base segment, this value is positive when 
pressure is directed downward and negative when 
directed upward. 


IP=2, value of PV is disregarded and a dead load 
analysis is carried out for the unit weights 
specified on the SEGMENT PROPERTIES cards. 


IP=3, PV is the magnitude of the uniformly 
distributed prestress pressure on the 
midsurface. Same sign convention as IP=1. 


IP=4, PV is the uniform temperature change in 
degree Celsius or Fahrenheit, depending on the 
units of ALPHA. (positive if the temperature 
rises above the reference temperature) 


IP=5, PV is the gradient of temperature across 
section in degrees per unit of thickness. 
(positive if the temperature rises above the 
reference temperature) 


IP=6, PV is the magnitude of uniform pressure 
distributed over horizontal’ projection, or 
snow load, in psf 
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If IP=7, PV is the magnitude of the liquid weight 
density in pcf. 


WHT = height of liquid above the vertex of a cone or 
height of liquid above the inner edge of a sphere in 
feet. Value is ignored for load types other than 
liquid pressure. 


PSF(1) = magnitude of externally applied horizontal 
force at the top of the segment, lbs/ft. 


PSF(2) = magnitude of externally applied moment at the 
top of the seoment, Et-lbs/ALe. 


PSF(3) = magnitude of externally applied horizontal 
force at the bottom of the segment, lbs/ft. 


PSF(4) = magnitude of externally applied moment at the 
bottom of the segment, ft-lbs/ft. 


PSF(5) = magnitude of externally applied vertical force 
at the top of the segment, lbs/ft. 


PSF(6) = magnitude of externally applied vertical force 
at the bottom of the segment, lbs/ft. 


NOTE?: 

The PSF forces are forces and moments which, if 
necessary, are to be applied IN ADDITION TO the 
Gistributed loading effects identified by the PV values. 

Prestressing effects are generally simulated as 
distributed loads but cable anchorages give rise to 
concentrated loads which are treated as PSF forces. 
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FLEXSHELL Input and Output Files 


for the Intze Tank Problem 
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FLEXSHELL Program Listing 


CO nae | es, 


bl Pe 


- 
OWMIMHUN EWN = 


~— 
N= 


AAANKHDAANANAAAANANAANANANANANANANANHANANNANAANAANANAANANANANAANANANANANANAANANNHAANAANANHNANANAAANAANNH 


aes FLEXSHELL ¥Bz 


THIS PROGRAM WAS ORIGINALLY DEVELOPED BY D.W. MURRAY 
AND A.M. ROHARDT IN DECEMBER, 1976 (REVISED IN MAY, 1977) 
FOR THE FLEXIBILITY ANALYSIS OF SEGMENTED AXISYMMETRIC 
SHELLS SUCH AS THE ‘SHORT’ CYLINDER, ‘LONG’ SPHERE, AND 
BASE SEGMENTS, SUBJECTED TO DEAD LOAD, UNIFORM PRESSURE 
GR PRESTRESSING, AND THERMAL STRESSES. REDUNDANT FORCES 
ARE APPLIED TO THE INDIVIDUAL SEGMENTS TO ESTABLISH THE 
REQUIRED GEOMETRIC COMPATIBILITY. 


THE CAPABILITY OF THE PROGRAM TO HANDLE &A WIDER 
VARIETY OF PROBLEMS WAS INCREASED IN SPRING, 1983 BY 
N. HERNANDEZ. IN ADDITION TO THE EXISTING SEGMENTS, THE 
“SHORT’ SPHERE, CONE, INVERTED SPHERE, AND INVERTED CONE 
WERE ADDED. THE ‘LONG’ SPHERE WAS REPLACED. TwO LOAD CASES 
WERE ALSO INCLUDED: THE SNOW LOAD WHICH IS A UNIFORM 
PRESSORE OVER THE HORIZONTAL PROJECTION OF THE SEGMENT, 
AND THE LIQUID PRESSURE LOADING. CONSEQUENTLY, THE 
FOLLOWING OPERATIONS WERE COMPLETELY RECODED: 


CALCULATION OF THE MEMBRANE STRESSES. 
CALCULATION OF THE EFFECTS OF A VERTICAL EDGE LOAD. 
CONSTRUCTION OF THE FLEXIBILITY MATRIX. 
CALCULATION OF THE MEMBRANE DISPLACEMENTS. 
CALCULATION OF THE FINAL STRESS RESULTANTS AND 

* OISPLACEMENTS. 


WeUun— 


sss NOTATION #4 


IT SEGMENT TYPE 

1 CYLINDER 

2 SPHERICAL DOME 

3 CONICAL DOME 

4 BASE ON ELASTIC FOUNDATION 
5 INVERTED SPHERE 

6 INVERTED CONE 
IP TYPE OF LOADING 

1 INTERNAL PRESSURE 

2 DEAD LOAD 

3 PRESTRESS 

4 UNIFORM THERMAL STRAIN 

s GRADIENT THERMAL STRAIN 

6 UNIFORM LOAD OVER 4&4 HORIZONTAL PROJECTION 


LIQUID PRESSURE 


GEOMETRIC VARIABLES 


SEG. TYPE CYLINDER SPHERE CONE BASE 
T THICKNESS THICKNESS THICKNESS THICKNESS 
R RADIUS RADIUS SEMI-VERTEX ANGLE BASE STIFFN 
H LENGTH OUTER ANGLE OIST. FR. VERTEX OUTER RADIUS 


TO LARGE END 


HO bald INNER ANGLE OIST. FR. VERTEX INNER RADIUS 


TO SMALL END 


INDECES 
NF=NUMBER OF SEGMENT EDGE FORCES 
NSEG=NUMBER OF SEGMENTS 


MAIN ARRAYS 
IR(*,1)3REDUNDANT FLAG TOP OF ELEMENT 
IR(*,2)SREOUNDANT FLAG BOTTOM OF ELEMENT 
IDFSIDENTITY OF UNKNOWN FORCES AT TOP AND BOTTOM OF ELEMENTS 
PBF=PARTICULAR SOLUTION BASE FORCES 
PSD=PARTICULAR SOLUTION EDGE DISPLACEMENTS 
PARD=PARTICULAR SOLUTION INCOMPATIBLE DISPLACEMENTS 


PSFEPARTICULAR SOLUTION FORCES WHICH PRODUCE ADDITIONAL INCOMPAT 


-IBLE DISPLACEMENTS 


A=MATRIX ESTABLISHING GEOMETRIC COMPATIBILITY BETWEEN DEGREES OF 


FREEOOM 
EXTERNAL FUNCTIONS AND SUBROUTINES FOLLOW THE MAIN PROGRAM 
IN THE FOLLOWING ORDER: 


FUNCTIONS: SUBROUTINES: 

a PNA 4. PCYLIN 9. BASE 14. PCONE 

27 EN2 Ss. CYUIN 10. BSHAPE 15. CONE 

3. FN3 6. POOME 11. JINVER 16. MMKEL2 

a. FNS 7. DOME 12. SOL Lecdacte wie OL Th 
& PBASE 13. PFOR 18. ROWEX 


IMPLICIT REAL*#8(A-H,O0O-2) 

DIMENSION 1T(20),R(20),H(20) ,HO(20!},E(20),PR(20) ,ALPHA(20) ,PV(20}, 
* $(6,6),PSD0(6),F(80,80),TT(80,80),PARD(80O) ,PART(80), PBF(6,20), 

x FO(80O),SF(6),CVEC(4) ,XH(101),A(8C,80),TS(4,4) ,PS©&(20,6),BB‘(4,4), 
= RM1(101),RM2(101),RN1(101),RN2(101),EC(20,2) ,UW(20),TITLE(10) 
DIMENSION I7T(20),1R(20,2),1P(20;,10C0(18,2),10F(19,6€),ND1V¥(20) 

* ,I1BASE(6),I1VECT(4),SR(10),M1(100) ,M2(100),V(100),HR‘(101),W( 100) 
x ,WHT (20) 

DATA PI/3.1415926536/ ,RAD/S7.295779513/ ,IBASE/1,2,5,3,4,6/ 


READ(S5,1001) TITLE 

WRITE(6,2001) TITLE 

READ(5,1000) NSEG, IPRINT 

WRITE(6,2000) NSEG,IPRINT 

IF(NSEG.GT.20) Go To 399 

READS) 1000) SC1ALUOD), IRIT, 1), TROL. 2) NOIVGL EC Ch ABCC 20, 
x I=i1,NSEG) : 

WRITECE, 2100) (CLL ITCL) IRI, 6). IROL, 29 NOIV CL) CCCI) Eett,2), 
* 12=1,NSEG) 

NSEGI=NSEG~-1 


READ(5,1200) ((I1DCO(I,J),J21,2),1=1,NSEG1) 
WRITE(6,2200) ((IDCO(1,J),J=1,2),1=1,NSEG1) 

READ(S,1300) (1,T(1),R¢1),H(1),HO(1),E(1),PR(T),ALPHA(I),UW(I), 
* [21,NSEG) 

WRITE(6,2300) (1,T(1),R(1),H(I),HO(I),E(1),PR(1),ALPHA(I),UW(I1), 
s 1=1,NSEG) 

READ(S,1400) (1,I1P(1),PV(I),WHT(I),(PSF(I,J),J=1,6),1=1,NSEG) 
WRITE(6,2400) (1,1P(1),PVi1),WHT(1),(PSF(I,J)},H=1,6),1=1,NSEG? 
pa=PI/4 

R2=2.0 


RZ2=DSORT(R2> 


NR2NUMBER OF REDUNDANTS 
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C IDENTIFY DEGREES OF FREEDOM AND FORM A-MATRIX 


Crore er eee eee ee c meer ec eee cee eee ee eee eee ewe eee ere ce eer ecco ec ee cee ee ec ee 


c FORM IDF ARRAY 
00 30 J=1,6 

00 30 121,NSEG 
30 IDF(I,J)=0 


KVz0 
KOUNTEZO 
NRHZO 
DO SO I#1,NSEG 
DG 560 J31,2 
J22¢J-1 
IF(IR(I,J).E9.0) GO TO S90 
IDF(1I,J2)sKOUNT+1 
IF(IR(I,J).GT.0O) GOTO 40 
KOUNT=KOUNT#+1 
NRH=- i 
GOTO 41 

40 IDF(1,J2+1)=KOUNT*2 
KOUNT=KOUNT4+2 

4i IF(IT(1I).NE.3) GO TO 50 
IF(J.EQ.1) J225 
IF(J.EQ.2) J226 
KOUNTSKOUNT+i 
IDF(1,J2)=KOUNT 
KVSKV+1 

50 CONTINUE 

NF=KOUNT 
WR=2*NSEG1+KV/2+NRH 
IF(NF.GT.80) GO TO $99 


IF(IPRINT.EQ.0) GOTO 45 
WRITE(6,2435) 
00 35 1=1,NSEG 
WRITE(6,2440) (IDF(I1,J), Je1,6) 
35 CONTINUE 
2440 FORMAT(20I15) 
2435 FORMAT (///* THE IOF MATRIX IS ‘/) 
4s CONTINUE 


c FORM A MATRIX 
1021 
BDO 180 I=1,NR 
DO 180 J=1,NF 

180 A(I,J)=0.0 
BDO 300 I=1,NSEG1 
K=IDCO(1,1) 
L=I1o0C0(1,2) 

250 JISIDE(K,3) 
J2=1IOF(L,1) 
A(I0,JjJ1)=1. 
A(I0,JjJ2)=-1. 
IF(IDF(K,4).EQ9.0.0R.IDF(L,2).EQ.90) GOTO 260 
A(TO+1,JSit1)=1. 
A(IO#1,JS241)2-1. 


260 IF(IT(L).EQ.3) A(I10,J241)S4EC(L,1) 


IFCIT(K).EQ.3.AND.IT(L).EQ.3) AC(IO,J241)=-EC(K,2) 


1o=10+2 
IF(IDF(K,4).EQ.0.0R.I1DF(L, 
IFCIT( I). NE.3.OR.IT(L) .NE. 
JI=IDF(K,6) 

J2SIDF(L,5) 

A(I0,J1)=1 

A(10,J2)s- 

Io=!10+1 

300 CONTINUE 


).EQ0.0) 10=10-1 
) 


2 
3) GOTO 300 


IF(IPRINT.EQ.0) GO TO 351 

WRITE(6,2450) 

DO 350 I=1,NR 

WRITE(6,2500) (A(I,J),S=1,NF) 
3s50 CONTINUE 
2450 FORMAT (///* THE A CONNECTIVITY MATRIX IS%) 
2500 FORMAT(20F4.1) 


Cocrrcersercorzres ecrcerercec ecco e reer esr ene eee eee eee eee eee sere eee reese 
351 DO 355 N=:,NSEG 
DO 355 J=1,6 
385 PBF(J,N)=0.0 
OC 356 I=1,NSEG 
356 PBF(5,3])=PSF(1,5) 
c 
DO 395 N=1,NSEG 
c 
I0ov = 0 
CALL PFOR(IT(N),T(N),R(N),H(N),HO(N) ,WHT(N),E(N),PRON), UWC(N) 
= ,ALPHA(N),IP(N),PV(N),PBF(5,N),PBF,IDV,N) 
c 
ITN = IT(N) 
GOTO (360,361,362,363,364,365),ITN 
c E 
c CYLINOER 
c 
360 RN = RIN) 
RO = 1 
Go TO 370 
c 
c SPHERE 
- 
361 RN = R(N)*#OSIN(H(N)/RAD) 
RO = DSIN(HO(N)/RAD)/OSIN(H(N)/RAD) 
234850. 
IF(HO(N) .NE.O. ) Z1 = DCOS(HO(N)/RAD)/OSIN(HO(N)/RAD) 
22 = DCOS(H(N)/RAD)/DSIN(H(N)/RAD) 
Go TO 369 
c 
& BASE ON ELASTIC FOUNDATION 
x 
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362 


391 


ano 


365 


385 
395 


IF(HO(N).LT.1.0E-03) GO TO 391 
PSF(N,4)=SPSF(N,4)+PB8F(4,N) 
PSF(N,2)=PSF(N,2)+PBF(2,N) 
PSF(N,S)=PBF(S,N) 

Go TO 395 


CONE 


RN 
RO 
z1 
22 
Go TO 369 


H(N)*OSIN(R(N)/RAD) 
HO(N)/H(N) 
DTAN(R(N)/RAD) 


INVERTED SPHERE 


RN = R(N)*DSIN(HO(N) /RAD) 

Zi = -DCOS(H(N)/RAD)/DOSIN(H(N)/RAD) 
RO = 0. 

z72 2 0 


IF (HO(N).EQ.0.) GO TO 368 

RO = DSIN(H(N)/RAD)/DSIN(HO(N)/RAD) 
22 = -DCOS(HO(N)/RAD)/DSIN(HO(N)/RAD) 
GO TO 3693 


INVERTED CONE 


RN 
z1 
pS 3 Z1 

RO °. 

IF (HO(N).EQ.0.) GO TO 369 

RO = H(N)/HO(N) 
PSF(N,1)2PSF(N,1)+PBF(S,N)*#Z1 
PSF(N,2)=-PBF(S,N)#EC(N,1)*PBF(2,N)+PSF(N,2) 


HO(N)*OSIN(R(N) /RAD) 
-DTAN(R(N)/RAD) 


PSF(N,4)©SPSF(N,4)¢(PBF(5,N)*RO+PBF(G6,N))#EC(N,2)+PBF(4,N) 


DO 385 T=1,NSEG1 
IF(IOCO(1I,1).NE.N) GO TO 385 
L=1oco(I,2) 

RL = R(L) 


RL*DSIN(HO(L)/RAD) 
H(L) 
HO(L)*DSIN(RL/RAD) 


IF(IT(L).£O.2) RL 
IF(IT(L).EO.3) RL 
IF(CIT(L).EO.4) RL 
TFCIT(L). 80.58) RL RL#O0SIN(H(L)/RAD) 
IF(IT(L) .EQ.6) RL H(L)*DSIN(RL/RAD) 
IF(RL.LT.1.00-06) GO TO $998 


une ne 


IF(ITN.NE.1) PSF(L,1)SPSF(L,1)4(PBF(3,N)-PBF(S,N)*RO#Z2)=RN/RL 


PBF(S,L)=PBF(S,L)+(PBF(5,N)#RO+PBF(6,N))*RN/RL 
IF(IT(L).E9O.3) PSF(L,2)SPSF(L,2)-PBF(3,N)*EC(L,1) 
CONTINUE 

CONTINUE 

LOVe = rat 

IF(IPRINT.EO.0)GOTO 398 

WRITE(6,2550) 

DO 396 I=1,NSEG 

WRITE(6.2551) (PBF(J,1),J21,6) 

WRITE(6, 2552) 

DO 397 1I2=1,NSEG 

WRITE(6,2553) (PSF(I,J),J=1,6) 

FORMAT(///,’ ase PBF zezx ‘/) 

FORMAT(6E13.4) 

FORMAT(///,’ ee PSF xe 2) 
FORMAT(6E13.4) 


c CONSTRUCT AND ASSEMBLE ELEMENT FLEXIBILITY MATRICES 


14/ 


C AND INITIAL DISPLACEMENT VECTOR 
Corre err cre cere score ee ees orc eeee micieletiel= 
338 DO 400 1=1,NF 
PARD(1)=0.0 
DO 400 J=1,NF 
400 F(T, Jd) s0°.0 
¢€ 
DO 500 N=1,NSEG 
ITN = IT(N) 
IFLAG=0 
GOTO (401,405,410,470,4805,470),ITN 
401 CALL CYLIN(T(N),R(N),H(N),HO(N),E(N),PR(N),UW(N),S,TS,0,BETA, 
= IFLAG) 
PBT=PBF(5,N) 
CALL PCYLIN(T(N),R(N),H(N),HO(N),E(N),PR(N),UW(N),ALPHA(N),S,PSD, 
* JIP(N),PV(N),N,PSF,PBT) 
c 
IF(IPRINT.EQ.0) GO TO 420 
WRITE(6,260C) N,((S(1I,J),JS21,4)3,12=1,4! 
2600 FORMAT(///’ FLEXIBILITY MATRIX FOR CYLINDRICAL SEGMENT’, 
= (4E16.5)) 
c 
GOTO 420 
(3 
C SPHERICAL SEGMENT 
405 CALL DOME(IT(N),T(N),R(N),H(N),HO(N),E(N),PR(N),UW(N),S,TS 
x» ,ANG,ANGO,RLAM,IFLAG) ; case : 
PBT=PBF(5,N) 
CALL PDOME(ITI(N),T(N),R(N),H(N),HO(N),E(N)., PRON) ,UW(N) ,ALPHAI(N), 
= S,PSD,IP(N),PV(N),I1DV,PSF,ANG,ANGO,WHT(N),PBT,IFLAG,N) 
c 
IF (IPRINT.EQ.0) GO TO 420 
WRITE(6,2700) N,((S(1I,J),8=1,4),121,4) 
2700 FORMAT (///* FLEXIBILITY MATRIX FOR DOME SEGMENT’ ,14/ 
= (4E16.5)) 
GOTO 420 
c 
c ELASTIC FOUNDATION SEGMENTS 
410 CALL BASE(IFLAG,T(N),R(N),H(N),HO(N),E(N),PR(N) ,UW(N),BE,S,TS,D) 
CALL PBASE(T(N),R(N),H(N),HO(N),E(N),PR(N),ALPHA(N) ,UW(N),S,PSD, 
x IP(N),PV(N),N,PSF, PBF) 
c 
fe: SPECIAL ASSEMBLY FOR BASE ELEMENTS 
c 


DO 415 1=1,6 
ETRE VN Ta 
IF(L.EQ.0) GO TO 415 


ELZ 
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PARD(L)=PSD(I) 


412 J2#1,6 


K=IDFI(N,J) 


IF 


(K.EQ.0) GO TO 412 


F(L,K)2S(I,J) 
412 CONTINUE 
415 CONTINUE 


KE 


CIPRINT.EQ.0) GO TO SOO 


WRITE(6,2720) N,S 
2720 FORMAT(///’ FLEXIBILITY MATRIX FOR BASE SEGMENT’ ,14/ 


(6E16.5)) 
To 500 


C CONICAL SEGMENT 
470 CALL CONE(IT(N),T(N),R(N),H(N),HO(N),E(N),PR(N) ,UW(N),S,TS 
*,ANG,XM,RLAM, IFLAG) 
PBT=PBF(S,N) 
CALL PCONE(IT(N),T(N),R(N),H(N),HO(N),E(N),PR(N),UW(N),ALPHA(N), 


aa 


S$,PSO,IP(N),PV(N),IDV,PSF,ANG,WHT(N),PBT,IFLAG,N) 


(IPRINT.EQ.0) GO TO 420 


WRITE(6,2750) N,((S(I,J),321,4),121,4) 
2750 FORMAT(///°* FLEXIBILITY MATRIX FOR CONE SEGMENT’ ,14/ 


= 
c 


C ASSEMBLY OF FLEXIBILITY MATRIX AND DISPLACEMENT VECTOR 


(4E16.5)) 


420 DO 460 131,4 
L=IDF(N,1) 
IF(L.£€0.0) GO TO 4690 
PARD(L)=PSO0(1I) 

BDO 450 J=1,4 
K=IOF(N,J) 
IF(K.EQ9.0) GO TO 450 


F(L,K)=S(1,0J) 
450 CONTINUE 
460 CONTINUE 
500 CONTINUE 


ae 


(IPRINT.EO.0O) GO TO S08 


WRITE(6,2800) 
2800 FORMAT(///‘ ELEMENT FLEXIBILITIES AFTER ASSEMBLY *) 


bo 


sos 1=1,NF 


WRITE(6,2850) (F(I,J5),JS=1,NF) 
505 CONTINUE 
2850 FORMAT(6E12.4) 


c 


(PARD) 


c CONDENSE TO REDUNDANT FLEXIBILITY MATRIX AND DISPLACEMENT VECTOR 


508 bo 
oo 


520 I=1,NF 
S20 J=1,NR 


c=0.0 


510 K=1,NF 


510 C=C+F(1,K)¥A(J,K) 


520 wu 


(1,J)=2¢ 


650 1=1,NR 


c=0.0 


do 


610 K=1,NF 


610 C=C+A(1,K)*PARD(K) 
PART(I)=-C 


c 
DO 630 J=1,NR 
c=0.0 
DO 620 K=1,NF 
620 CzC+A(1,K)*FTT(K, J) 
630 F(I,J)=C 
650 CONTINUE 
c 
IF (IPRINT.EO.0O) GO TO 670 
WRITE(6,2900) 
2900 FORMAT(///‘ CONDENSED FLEXIBILITY MATRIX’) 
DO 660 I1=1,NR 
WRITE(E6€,2850) (F(1I,J),JS=1,NR) 
660 CONTINUE 
WRITE(6,3001) (1,PART(1),1=1,NR) 
3001 FORMAT(///‘ INCOMPATIBLE DISPLACEMENTS’/’ IDISP VALUE’ / 
z GUS ZEv3) 50.) 
Coreen cansverece ecrceee eer eee eee wee eee eee ee eee eee es 
c SOLVE FOR REDUNDANTS AND FIND SEGMENT END FORCES 
Ew ess eecerec ese em se id eect ee ad err ee ee eee eee eee e-ee 
670 CALL SOL(F,PART,80O,NR) 
c 
c FIND SEGMENT FORCES 
do 700 =1,NF 
c=0.0 
Do 690 J=i,NR 
690 CHC+A( J, 1)*PARTi J) 
700 FO(I)=sc 
c 
c WRITE TOTAL VECTOR OF SEGMENT END FORCES 
WRITE(6,3100? 
3100 FORMAT(///‘FORCES ON ENDS OF SEGMENTS’/’ SEG J ea 6x, 
= “FORCE” ) 
00 706 N=1,NSEG 
DO 70S J=1,€ 
L=IOFI(N,J) 
IF(L.EQ9.0) GO TO 705 
WRITE(6,3200) N,J,L,FO(L) 
705 CONTINUE 


706 CONTINUE 
3200 FORMAT(314,E13.5) 


c ceoereoecesneevenecseececeneeeree2e2 e220 e022 2222°2° 
c FIND AND OUTPUT SEGMENT STRESS RESULTANTS 
C------ Bia) a aba aiela) aibiete leis sta a aia) s! ainiw ewe ise iw. e\ele ae ome wie lehei~ a Sie ees (6 allots aie) 6 
c 

DO 900 N=1,NSEG 

IFLAG=1 

ITN = ITI(N?> 

IPN=IP(N) 

IF(CIT(N).EQ.3) GO TO 800 
c 
c FORM SEGMENT ENO FORCE CVECTOR 


IF 


(NDIV(N).GT.100) Go Tc 939 


DH=(H(N)-HO(N) )/DFLOAT(NDIV(N)) 
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NOIVIZNDIV(N) +1 
00 710 121,4 
SF(I)=PSF(N,1) 
L=IDF(N,I) 
IF(L.EO.0) GO TO 710 
SF(I)=FO(L)+*SF(I) 
710 CONTINUE 
SF(S)=PBF(5,N) 
SF(6)=PBF(6,N) : 


C WRITE INDIVIDUAL SEGMENT END FORCES 
WRITE(6,3300) N,(SF(I),121,6) 

da FORMAT (///" INDIVIDUAL END FORCES FOR SEGMENT ’,16,/(4E13.5)) 
GOTO (708,750,800,870,750,870),1TN 


Creer ce eww ew ecw cee we ee ee eee ee wee mew eee meee ec ce cece eee eee ee 


c STRESS RESULTANTS FOR CYLINDRICAL SEGMENTS 
Ce ces Sic eb wie ica wieich ssac Vesa el G)e Mal a eels ala cic aldlaiee ales s'sld ele sidisleieie da aitelea dae «aac 
70S CALL CYLIN(T(N),R(N),H(N),HO(N),E(N),PR(N),UW(N),S,TS,0,BETA, 
* IFLAG) 
ON1=0.0 
CM1=0.0 
CM2=2=0.0 
CNi=-PBF(S,N) . 
WP=SCN1#®R(N)*PR(N)/(E(N)D* TON) ) 
CN2=0.0 
WWHO.O 
WL=0.0 
RNSE(N)*T(N)/R(N) 
GOTO (711,712,711,713,714,714,715),IPN 
711 CN22-PV(N)*R(N) 
WPSPV(N)#®R(N)**#2/(E(N)*®T(N) ) #WP 
GoTo 721 
712 DNi=-T(N) *UW(N) *DH 
WW=EON1*R(N)*PR(N)/(E(N)*T(N)) 
GOTO 721 
713 WP=-R(N)*ALPHA(N)*®PV(N) WP 
GOTO 721 
714 CM1=(1.+PR(N))*®O2ALPHA(N) *PV(N) 
CM2=CMi 
Go TO 721 
715 WL=PV(N)®R(N)#*2/(E(N)*T(N)) 
721 DZ. *O0sBETA*+2 
00 730 1=1,4 
c=0.0 
BO 720 J=1,4 
720 C=C+TS(1,J)*SF(J) 
730 CvEC(I)=sc 


x20.0 

DO 740 t=1,NDIVI 
BXSBETA#X 
Bpc=o0co0sS(Bx) 
DS=DSIN(BX) 


W(L) SWP+DEXP(BX)*(CVEC(1)*0C+CVEC(2)*D0S)+DEXP(-BX)*# 
* (CVEC(3)=O0C+CVEC(4)*0S) +WWeDFLOAT(L-1)+WL*X 


RN1(L)SCN1+*DNI*DFLOAT(L-1) 
RN2(L)SCN2-(DEXP(BX)2*(CVEC(1)*0C*#CVEC(2)2*0S)+ 
=DEXP(-BxX)*(CVEC(3)*DC+CVEC(4)*DS)+WL*X)#RN 
RM1(L)SDEXP(BX)*(D*CVEC(2)*0C-D*CVEC(1)*D0S)+DEXP(-BxX)*(D= 
* CVEC(3)*0S-D*CVEC(4)*0C) : 
RM2(L)=PR(N)#RM1(L)+CM2 
RM1(L)=RM1(L)+CM1 
KH(LY=X 

740 X=X+DH 


WRITE(6,4000) N,(L,XH(L),RNI(L),RNZ2(L),RM1(L),RM2(L),LE1,NDIV1) 
WRITE(6,4005) (L,XH(L),W(L),L=1,NDIV1) 
GOTO 900 


Giaies <a ae aa a sta ate le © (a=) Sielerarsisleseis e's ere eee et eee ed 


750 CALL DOME(IT(N),T(N),R(N),H(N),HO(N),E(N),PR(N) ,UW(N),S,TS 
, ANG, ANGO,RLAM,IFLAG) 


” 


x=O. 

OX=DH/RAD 

D=ZE(N)2T(N)**3/(12.2(1.°PR(N)*®*2)) 
c 

DO 767 1=1,4 

C=0.0 

00 766 J=1,4 
766 CHeC+TS(1,J)*SF(Ji 
767 Cveci1)=C 
c 

00 780 L=1,NOIVi 
c 

IFCITIN) .EQ.2) PHI = X+ANGO 

NE PWC NDE ./S)) PHI] = ANG-xX 

CAOX=0COS(PHI) 

SAQOX=DSIN(PHI) 
c 

oco = DCOS(RLAM*PHI } 

OS1 = DSIN(RLAM*PHI) 

EP = DEXP(RLAM=PHI) 

EM = DEXP(-RLAM*PHI) 

TH1 = EP*(DCO*DSI) 

TH2 = EP*(DCO-OSI) 

TH3 = EM*(DCO+DS!1) 

TH4 = EM*(DCO-OSI) 

PHI1= EP*OCO 

PHI2= EP#DSI 

PHI3= EM*DCO 

PHIG= EM=DSI 
c 

WP=0 

RM1(L) = O. 

RM2(L) = 0 

RN2(L) = FN2(PV(N),SF(S5S) ,UW(N),R(N),T(NE,IP(ND, 

x IT(N),PHI,A8NGO,ANG,WHT(N), IDV) 


GOTO (753,782,753,754,760,753,753),I1PN 
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753 RNI(L) = FN1(PV(N),SF(S),UW(N),R(N),T(N),IP(N), 
x IT(N),PHI,ANGO,ANG,WHT(N),IDV) 
WP = -R(N)*SAOX*(RN2(L)-PR(N)*RNI(L))/(E(N)*T(N) ) 
GO TO 755 
754 WP=-ALPHA(N)*R(N)*SAOX®PV(N) 
GO To 755 
760 RMi(L) = PV(N)*ALPHA(N)*D2(1.+PR(N))/T(N) 
RM2(L) = RM1I(L) 
c 
755 wW(L) = WP + R(N)*SAOX*RLAM*®(CVEC(1)*TH2 + CVEC(2) 
A * *TH1 - CVEC(3)*TH3 + CVEC(4)*THS)/(E(N)#*T(N)) 
RN2(L) = RNZ(L) - RLAM#*(CVEC(1)#TH2 + 
E * CVEC(2)*TH1 - CVEC(3)*TH3 + CVEC(4)*TH4) 
758 IF(PHI.EQ.0.) GO TO 761 
RN1T(L) = RNI(L) = (C&OX/SAOX)*(CVEC(1)*PHI1 
‘ * + CVEC(2)*PHI2 + CVEC(3)*PHI3 *+CVEC(4)*PHIG) 
761 RM1(L) = RMI(L) + 0.5#R(N)*(-CVEC(1)#£TH1 + CVEC(2)#TH2 
* + CVEC(3)*TH4 + CVEC(4)*TH3)/RLAM 
c 
RM2(L) = RM2(L) + PR(N)*RM1(L) 
¢ e 
771 X=X+DX 
XH(L)=SDFLOAT(L-1)*0H 
780 CONTINUE 
c 
WRITE(6,4500) N,(L,XH(L),RN1(L),RN2(L),RM1(L),Q2M2(L),L=21,NDIV1) 
WRITE(6,4005) (L,XH(L),W(L),L=1,NDIV1) 
GOTO 300 
c 
| erreere 
C STRESS RESULTANTS FOR BASE-SLAB ELEMENT 
FREI ES TO SER ie CIE SO pL a te 


c FORM SEGMENT END FORCE VECTOR 


800 IF (NDIV(N).GT.100) GO TO 399 
NDIViSNDIV(N) +1 
DO gsoOS 1=1,6 
SF(I)=PSF(N,I1) 
L=IDF(N,I1) 
IF(L.EQ.0) GOTO 805 
SF(I)SFO(L)+SF(I) 
805s CONTINUE 
c 
c WRITE INDIVIDUAL SEGMENT END FORCES 
WRITE(6,5000) N,(SF(I),1=1,6) 
5000 FORMAT(///’ INDIVIDUAL END FORCES FOR’, 
* * SEGMENT’,16,/(4E13.58)) 
c 


c COMPUTE STRESS RESULTANTS 


815 C=xC+TS(1,J5)*SF( JJ) 
820 CVEC(I)=C 
HT=H(N) 
(= 
IFLAG=2 
DO 825 11,4 
DO 825 J=1,4 
825 BB(1,J)=0.0 
DO 860 L=1,NDIV1 
CALL BASE(IFLAG,T(N),R(N),HT,HO(N),E(N),PR(N),UW(N),BB,S,TS,D) 
DO 830 1=1,4 
cc=0.0 
DO 835 J=1,4 
835 CC=CC+BB(1,J)*CVEC( J) 
830 SR(I)=sCC 
HR(L)=HT 
DH=(H(N)-HO(N))/DFLOAT(NDIV(N) ) 
HT=HT-DH 
c . 
RM1(L)=SR(2) 
RM2(L)=SR(1) 
V(L)=SR(3) 
W(L)=SR(4) 
GOTO (840,841,860, 860,860,860,840),IPN 
840 W(L)EW(L)*#PVOIN)/RON) 
GOTO 860 
841 W(L)SW(L) UWE ND *ETIN) /RON) 
GOTO 860 
860 CONTINUE 
WRITE(6,5500) N,(L,HR(L!,RM1(L),RM2(L),V(L),Le1,NDIV1) 
WRITE(6,5505) (L,HR(L),W(L),L=1,NDIV1) 
GO TO 900 
c 
Cerrtcree eer eee eer eee ee eee ad coer eee er we we ee oe eee eH ee ee BO 
c STRESS RESULTANTS FOR CONE SEGMENTS 
Cicieiaje = Sc eselesieceeies + ame es etc cre ee cece ee 
870 CALL CONE(IT(N),T(N),R(N),H(N),HOCN),E(N),PR(N),UW(N),S,TS 
*,ANG,XM,RLAM, IFLAG) 
c 
po 872 1=1,4 
¢ f= 30-. 
po 871 J=1,4 
871) SCE. TS: CHASE CJ): 
a72 GVEG(T)) = 7¢ 
c 
X=O0. 
D = E(N)*#T(N)##3/(12.2¢1.-PR(N)*®*2)) 
i > 
po 890 L=1,NDIV1 
c 


IFLAG=1 


CALL BASE(IFLAG,T(N),RON),H(N),HO(N),E(N),PR(N),UW(N),BB,S,TS,D) 


DATA IVECT/2,5,4,6/ 
DO 820 1=1,4 

cz0.0 

DO 815 J=1,4 
JJZIVECT( J) 


nunno 


ooo 
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soul he ts 


679 
680 
681 
682 
683 
684 
68s 
686 
687 
688 
689 
630 
691 
692 
693 
634 
695 
696 
697 
698 
698 


703 
704 


706 
7107 


c 


RM2(L) = 0. 


Y= xX + HO(N) 
IF(ITN.EQ.6) NSW CUR) = K 
IF(Y.EQ.0.) GO TO 884 


XI & 2.*®RLAM*Y*#*0.5 
CALL MMKEL2(XI,BER,BEI,XKER,XKEI, 
= OBER ,OBEI, DKER, DKE!I ) 


Gi = XI*DBER - 2.*PR(N)*BER 
G2 = XI#DBEI - 2.#PR(N)*BEI 
G3 = XI*#DKER - 2.*PR(N)*XKER 
G4 = MI*OKEI - 2.*PR(N)*XKEI 
Ui = Gi + 4.=PR(N)*BER 

U2 = G2 + 4.#PR(N)*BEI iS 
U3 = G3 + 4.*PR(N)*XKER 

U4 = G4 + 4. 2PR(N)*®XKEI 

Vie 2.*BER + PR(N)*XI*DBER 
V2 = 2.*BEI + PR(N)*XI*DBEI 
V3 = 2.®xXKER + PR(N)*XI*DKER 
V4 = 2.*XKEI + PR(N)*X120KEI 


880 GOTO (874,874,874,875,876,874,874),IPN 


874 


1001 

2001 
1000 
2000 


2100 
1200 
2200 
1300 
2300 


1400 
2400 


4000 


4005 


4500 


4510 


s500 


5505 


RNI(L) = FN3S(IP(N),IT(N),Y,ANG,H(N),HO(N), 


* WHT(N),PV(N),PBF(S,N),UW(N),T(N), IOV) 
RN2(L) = FNS(IPYN),IT(N),Y,ANG,H(N),HO(N), 
* WHT(N),.PV(N),PBF(S,N),UW(N),T(N)) 


WP = -Y*DSIN(ANG)*®(RN2(L)-PR(N)#RNI(L))/(E(N)#T(N)) 
GO TO 879 

WP = PV(N)#ALPHA(N)*Y*DSIN( ANG) *E(N)*T(N) 

GO TO 8793 

RM1(L) = PV(N)*ALPHA(N)*D*(1.+PR(N))/T(N) 

RM2(L) = RM1I(L) 


W(t) = (WP + 0.S#*DSIN(ANG)*®(CVEC(1)#G1 + CVEC(2)*G2 + 
= CVEC(3)*G3 + CVEC(4)*G4))/(E(N)*TIN) ) 


RNT(L) = RN1(L) - (CVEC(1)*BER + CVEC(2)*BEI + 
*CVEC(3)*XKER + CVEC(4)*XKEI)/Y 


RN2(L) = RNZ2(L) © O.5*XI#(CVEC(1)*0BER + CVEC(2) 
**DBEI + CVEC(3)*DKER + CVEC(4)*DKEI)/Y 


RM1(L) = RMI(L) - TC(N)*®(CVEC(1)*#U2 - CVEC(2)*U1 + 
#CVEC(3)*U4 - CVEC(4)#U3)/(2.*xKM2222Y) 

RM2(L) = RM2(L) - T(N)*®(CVEC(1)2V2 - CVEC(2)*V1 + 
=CVEC(3)*V4 - CVEC(4)#V3)/(2.*xXMB222Y) 

xX = Kk + DH 

XH(L) = DFLOAT(L-1)*0H 

CONTINUE 


WRITE(6,4510) N,(L,XH(L),RNI(L),RN2(L),RM1(L),RM2(L),L=1,NDIV1) 
WRITE(6,400S) (L,XH(L),W(L),L=1,NDIV1) 

CONTINUE 

STOP 

WRITE(6,3000) 

FORMAT(’ STOP FOR PROGRAM DIAGNOSED INPUT ERROR ‘*) 

stop 


FORMAT(10A8) 

FORMAT(‘1’,10A8//) 

FORMAT(514,2F10.4) 

FORMAT(‘1’,° sare OUTPUT FOR FLEXIBILITY ANALYSIS OF SEGMENTED’, 


= ' SHELL *#=8'/, 

= 0 NUMBER OF SEGMENTS S72 18/7 

* 4 IPRINT Ea a Whee 

* * SEG TYPE IR UR MONDIVG | 4X BOs 7h, CeCe) 
FORMAT(14,415,2F10.4) 

FORMAT (214) 

FORMAT(///,° CONNECTIVITY MATRIX’/(5X,214)) 
FORMAT(14,F6.0,F12.0,F8.0,5F10.0) 

FORMAT(///‘' GEOMETRIC PARAMETERS'/’ SEG’,4X,’THICK’,3X,’RADIUS' 


*,3X,’L OR ANG‘’,4X,‘ANGO’,7X,‘MODULUS’,6X,’°P RATIO’,2X, ’THERMCOEF’, 
* 3X, °WEIGHT’/ (16,F8.3,F12.3,2F10.3,613.4,F10.3,€13.4,F10.3)) 
FORMAT(214,8F10.0) 

FORMAT(///’ PARTICULAR SOLUTION INPUT INFORMATION’ / 

#7 ° SEG TYPE VALUE’ ,4X,’LIO HT’,2X,’TOP SHEAR’ ,4xX,’TOP MOM’, 

= 6X,’BOT SHEAR‘ ,4X,’BOT MOM’ ,4X,’TOP FORCE’ ,4xX,’VERT FORCE’/ 

eP ORCAS, F1OlS PIS 2S GETS ..5) 


FORMAT(///’ %*2 OUTPUT FOR CYLINDRICAL SEGMENT’ ,14,' == ='/ 
=’ POINT COORD N1°,12X,7N2°,12K,°M1°, 12K, °M2° ,/ 
2 (16,.F10.4,4E14.5)) ; 

FORMAT(///' *** HORIZONTAL DISPLACEMENT **2 ‘/ 
+ee) POINT COORD w ’/(16,F10.4,E14.5)) 

FORMAT(///’ = * QUTPUT FOR DOME SEGMENT’,14,° #==// 
=’ PCINT ANGLE Nags T2Kn N22 2K Mioe 12k5 OM 278,07 
2 (16,F10.4,4614.5)) 

FORMAT(///’ *** QUTPUT FOR CONE SEGMENT’,14,’ *22/'/ 
x’ POINT COORD Nive di2Ke ON2 AA 2k Mle kM 2 0 0/ 
* (16,F10.4,4E14.5)) 

FORMAT(///’ **= OUTPUT FOR BASE ELEMENT’,14,’ eax’/ 
<a POINT COORD Miiceutie kine M2 alia Xaver iit 


* /(16,F10.4,23814.8)) 

FORMAT(///’ *** VERTICAL DISPLACEMENT ***°*/ 
= ‘POINT coorD w‘/(16,F10.4,E14.5)) 
END 


CEASE SSAA SSAA KSLA TATA STATS A SRT TEARSTK SKA SLES AAT ASSEATTSHESE TE SLELEE 


FUNCTION FN1(PV,PBT,UW,R,T,IP,I1T,PHI,ANGO,ANG,WHT, IDV) 


c THIS FUNCTION IS USED FOR N1 DOME STRESS RESULTANTS 


5 


IMPLICIT REAL*8(A-H,0-2) 


FN1=0.0 

CATE TA 

GAMMA = ANGO 

PEARS NE SS) GOeT0!-s = 
CUPS Heise 


GAMMA = ANG = 
GOTO(10,,20, 10, 100; 100, 15,30), DP 
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c 


CSESSSSSSSESARAEKTRASSTSEHRSTTEFTETKAKTSRTEKASETEKT ESS ESETT STE TASES AKHKTEKLCSETEEETE 


c 


CESS ERE EKETAK EET A TST TAREE KKK STEAK ETE SATE STAAL ESHA SELES ETE EEE 


c 


CHEER RETEST REK SEER T TREAT T SKS SATS KTS SH ALERT ESSEAARAE SLE TARE TELE 


SUBROUTINE PFOR(IT,T,R,H,HO,WHT,E,PR,UW,ALPHA,IP,PV,PBT,PBF,IDV,N) 


c 
¢ 
c 
c 
c 


c 


10 FN1#-0.5*PV#R 
IF(ANGO.LE.1.06-03) RETURN 
FN1 = FN1#(1.-(DSIN(GAMMA)/DSIN(PHI))##2) 
GO TO 100 
15 FN1S-0.5*PVeR : 
IF(ANGO.LE.1.0€&-3) RETURN 
FN1 = C1*FN1*(1.-(DSIN( GAMMA) /DSIN(PHI))2##2) 
GO TO 100 
20 FN1 = -0O.S5*UWST#R 
IF (PHI.EQ.0.) RETURN 
FN1 © C1*2.2FN1*(DCOS (GAMMA) -DCOS(PHI))/DSIN(PHI)*##2 
GO TO 100 
30 «6oIF(PHI.EQ.0.) RETURN 
CONST1 = (DCOS(PHI)*#*#3-(D0COS(GAMMA)2#=#3))/DSIN(PHI)*#2 
CONST2 = (DSIN(GAMMA)/DSIN(PHI))#*2 
FN1 = -C1*PVER* (0.5% (WHT+C1*R*DCOS (GAMMA) )*(1.-CONST2) 
= +C1*R*CONST1/3.) 
100 IF(PHI.EQ.0.) RETURN 
IF(1IDV.EQ.0) RETURN 
FN1 = FN1 - PBT*OSIN( GAMMA) /DSIN(PHI)#*2 
RETURN 
END 


FUNCTION FN2(PV,PBT,UW,R,T,IP,IT,PHI,ANGO,ANG,WHT, IDV) 


THIS FUNCTION IS USED FOR N2 DOME STRESS RESULTANTS 
IMPLICIT REAL#8(A-H,0-2Z) 
FN220.0 
Crests 
GAMMA = ANGO 
TFECITEESC..2) Go To 5 
Gyte ates 
GAMMA = ANG 
S GOTO(10,20,10,100,100,15,30),IP 
10 FN25-0.S5*PVER 
IF(ANGO.LE.1.0E-03) RETURN 
FN2 = FN2*(1.+(DSIN(GAMMA)/OSIN(PHI))*#2#2) 
GO TO 100 
15 FNZ2 = -Cx.5#PVeR 
IF(PHI.EQ.0.) RETURN 
FN2 = C#0.S5*PV*R2(1.-(DSIN( GAMMA) /OSIN(PHI))2*2 
* - C*2.2DCOS(PHI)*22) 
GO TO 100 
20 FN2 = -C#0.5*UW*sT#ER 
IF (PHI.EQ.0.). RETURN 
FN2 = UW*tT*R*((DCOS(GAMMA)-DCOS(PHI))/ODSIN( PHI )*#*2 
* - C#DCOS(PHI)> 
Go To 1600 
30 IF(PHI.EQ.0.) RETURN 
CONST1 = (DCOS(PHI)#*3-DCOS(GAMMA)#*3)/DSIN(PHI)*#*2 
CONST2 = (OSIN(GAMMA)/DSIN(PHI))**#2 
FN2 = -C#PV2R2(0.5% (WHT+C#R*DCOS (GAMMA) )#(1.+CONST2) 
= -C#R*(CONST1/3.+DCOS(PHI))) 
100 IF(PHI.EQ.0.) RETURN 


IF (I1DV.EQ.0) RETURN 

FN2 = FN2 + PBT*DSIN(GAMMA) /DSIN(PHI)**2 
RETURN 

END 


FUNCTION FN3(IP,IT,Y,8NG,H,HO,WHT,PV,PBT,UW,T, IDV) 


THIS FUNCTION IS USED FOR Ni CONE STRESS RESULTANTS 
IMPLICIT REAL*8(A-H,0-Z) 
FN3 = 0. 
IF(Y.£Q.0.) RETURN 
Y1 = HO 
Cas. ae 
IF(IT.NE.6) GO TO § 
Give sic 
Yi =H 
S Chee -Cis0. S¥ Vee 2=Viite2 vy 
GOTO(10,20,10,100,100,10,40),1P 
10 FNS = PV*C#DTAN( ANG) 
GO TO 100 
20 FN3 = UW*TsC/DCOS(ANG) 
Go TO 100 
60 FN3 = -PV=Y*DTAN(ANG)*(3.*WHT#(1.-(Y1/Y)=*2) 
ra 4C0422.*Y*DCOS(ANG)*(1.-(Y1/Y)**3))/6. 
100 IF(IDV.EQ.0) RETURN 
FN3 = FN3 - PBT®Y1/(Y*DCOS( ANG) ) 
RETURN 
END 


FUNCTION FN4(IP,IT,Y,ANG,H,HO,WHT,PYV,PBT,UW,T) 


THIS FUNCTION IS USED FOR N2 CONE STRESS RESULTANTS 
IMPLICIT REAL#8(A-H,0-7Z) 
FNG = 0. 
Gat G7 ta 
Cc = -Y*OTAN(ANG) 


TROT TEORE eC Fa 
GOTO(10,20,10,100,100,30,40),IP 


10 FNS = -C#PV 
RETURN 3 
20 FNS = UWs=T*C#C1*DSIN( ANG) 
RETURN 
30 FNS = PV*¥C*¥C1#DSIN( ANG) 2#¥2 
RETURN 
40 FNG = -PV*Y*DTAN(ANG)*(WHT+C1*Y*OCOS( ANG) ) 
100 RETURN 
END 


THIS SUBROUTINE COMPUTES PARTICULAR SOLUTION EDGE FORCES (PBF) 
IMPLICIT REAL*8&(A-H,0-2Z) 
DIMENSION PBF(6,20) 


SELECT SEGMENT TYPE 
Go TO (50,500,600,700,500,700),IT 


CYLINDRICAL SEGMENTS 
50 GO TO ($00,100,900,900,650,900,300),IF 
DEAD LOAD 


100 PBF(6&,N)=PBF(6,N)+T*UW*H 


c 


RETURN 
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C SPHERICAL SEGMENTS 
500 ANG=H/57.295779513 
ANGO=HO/S7.235779513 
PHISANG 
Cae rarer. 
IFe(IT.&¢.2) GO TO 505 
PHI=ANGO 
Caan 
SOS RN1 = FNI(PV,PBT,UW,R,T,IP,IT,PHI,4NGO,ANG,WHT, IDV) 
GOTO (510,510,510,900,650,510,510) , IP 


510 PBF(6,N) = PBF(6,N) - RN1*OSIN(PHI) 
PBF(3,N) = PBF(3,M) + Ct#RN1*DCOS(PHI) 
RETURN 

c 

c BASE ON ELASTIC FOUNDATION 

600 GOTO (900,900,900,900,650,900,900) ,I 
c ° 


c THERMAL GRADIENT 

650 (idan Ties 
YECIT “GT. 3) Cz-1. 
PBF(4,N)2PBF(4,N)-C#E*PV®ALPHA#T#23/(12.%(1.-PR)) 
PBF(2,N)2-°-PBF(4,N) 
RETURN 

c 

C CONE SEGMENTS 

700 ANG2R/57.295779513 ° 
Yeu 
1352 jolie Sagi Is 

= DP OCET Ne, 6) GO TO 705 
Y = HO 
Ci 2 <1. 
7OS RNi = FN3S(IP,IT,Y,ANG,H,HO,WHT,PV,PBT,UW,T, IOV) 


GOTO (710,710,710,900,650,710,710),IP 

710 PBF(6,N) = PBF(6,N) - RN1*DCOS(ANG) * 
PBF(3,N) = PBF(3,N) + C1*RN1#O0SIN( ANG) 

$00 RETURN 
END 

CAHSEE RSK SKE RRAKSE RESTA TERRES TAA KSEE TERT ARETE TALE sLs esses 
SUBROUTINE CYLIN(T,R,H,HO,E,PR,UW,F,TT,0,BETA,IFLAG) 

c THIS SUBROUTINE COMPUTES THE CYLINDER FLEXIBILITY (F) 

c AND B MATRICES (TT). 
IMPLICIT REAL*8(A-H,O-2) 

4) 


DIMENSION F(6,6),TT(4, ,TA(4,4) - 
c 
D=EeT2*3/(12.%(1.-PR**2)) 
BETAS(3.4%(1.-PR#*2)/(R*T)#*#2)22.25 
PHI 1=0EXP(BETA*H) *0COS(BETA#H) 
PHI 2=DEXP(BETA*H)*DSIN(BETA#H) 
PHI3Z=DEXP(-BETA*H) *D0COS(BETA#H) 
PHIG=DEXP( -BETA*#H) ®0SIN(BETA#H) 
THI=PHI1+PHI2 
TH2Z=PHI1-PHI2 
TH3=PHI34PHI4 
TH4=PHI3-PHIG 
c 
TT(1,1) = +2. *D0*BETA##3 
TT(1,2) = 2.8D08BETA#3 . 
TT(1,3) = 2.*D*BETA®#3 
TT(1,4) = 2.*D*BETA#*3 
¢ 
Tira ciy Sy OF. 
TT(2,2) = +2.3D#BETA#2#2 
Jone ae =. Or 
TT(2,4) = 2.*DSRETA#2 
c 
TT(3,1) = TH1#2.*D0*BETA**3 
TT(3,2) =°TH2*2.*D*BETA*®#3 
TT(3,3) =-TH&a42. *O0*BETA*#3 
TT(3,4) =-TH34*2.*D*BETA#*3 
c 
TT(4,1) =-PHI24*2.*D*BETA®*2 
TT(4,2) = PH11*2.*D*#BETA®*2 
TT(4,3) = PHI4*2.*D*BETA=*2 
TT(4,4) =-PHI3*2.*D*BETA®=2 
CALL TTINV(TT,H? 
c 
IF(IFLAG.NE.O) RETURN 
c 
TAT fst 10 
TA(1,2)=0.0 
TA(1,3)=1.0 
TA(1,4)=0.0 
TA(2,1)=8ETA 
TA(2,2/=BETA 
TA(2,3)=-BETA 
TA(2,4)=BETA 
TA(3,1)=PHI1 
TA(3,2)=PHI2 
TA(3,3)=PHI3 
TA(3,4)=PHI4 
TA(4,1)=TH2*BETA 
TA(4,2)=THI*#BETS 
TA(4,3)=°TH3*BETA 
TA(4,4)=THO*BETA 
c 
DO 100 121,4 
DO 100 yv=1,4 
c=0.0 
pO 8O K=1,4 
Fe) C=C+TA(I,K)*#TT(K, J) 
100 F(1l,J)sc 
(= 
RETURN 
ENO 


CRAESTECESAATSRAAA CTSA SALSA CLASS TAR EASES TAR SE REA TS SETS EAA ELELERKRE LESS ATE 
SUBROUTINE DOME(IT,T,R,H,HO,E,PR,UW,F,TT,ANG,ANGO,RLAM,IFLAG) 

THIS SUBROUTINE COMPUTES FLEXIBILITY MATRICES (‘F+ FOR 

A COMPLETE OR TRUNCATED SPHERE. THE FLEXIBILITY MATRIX 

IS REDUCED TO A TWO BY TWO MATRIX FOR A COMPLETE SPHERE. 


nano 


IMPLICIT REAL*8&(4-H,0-2) 
DIMENSION F(6,£&),TT(4,4),TA(G,4) 


£23 


“ty 


oasst 


~ : ' he 
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“ear ‘oars +o enh 


Mab Md re i |, 1 + 5. : 


c 
ANG2H/57.295779513 
ANGO=HO/57.285779513 
RLAM=(3.4(1,.°PR¥*2)*(R/T)*¥*%2) 820.25 
PHI10 = DEXP(RLAM*ANGO) *DCOS(RLAM*ONGO) 
PHI20 = DEXP(RLAM*ANGO)*0SIN(RLAM*ONGO) 
PHI30 = DEXP(-RLAM*ANGO) *DCOS(RLAM*ANGO) 
PHI&0 = DEXP(-RLAM*ANGO)*DSIN(RLAM*ANGO) 
PHI11 = DEXP(RLAM#ANG)*DCOS(RLAM*ANG) 
PHI21 = DEXP(RLAM*ANG)*O0SIN(RLAMZANG) 
PHI31 = DEXP(-RLAM*ANG)*DOCOS(RLAM*# ANG) 
PHI41 = DEXP(-RLAM*ANG)*OSIN(RLAM*ANG) 
THI1O = PHI10 + PHI20 
TH20 = PHI10 - PHI20 
TH30 = PHI30 + PHI40 
TH40 = PHI30 - PHI40 
THT = PHI11 + PHI21 
TH21 =s PHI11 - PHI21 
TH31 = PHI31 + PHI41 
TH41 = PHI31 - PHI41 
IF (HO.NE.O.) Go TO 10 

c 


c INITIALIZE THE MATRICES TO ZERO. 


I=1,4 


DO 5 J=1,6 
TACT O): = Oh 


Ss TT(I,J) #& 0. 
GO TO 15 
10 TT(1,1) = PHI10/(-GSIN(ANGO) ) 
TT(1,2) = PHI20/(-DSIN(ANGO} ) 
TT(1,3) = PHI30/(-OSIN(ANGO) ) 
TT(1,4) = PHE@O/(-DSIN(ANGO) ) 
c 
TT(2,1) = -THIO#(-.5*R/RLAM) 
TT(2,2) = TH202(-.5#R/RLAM) 
TT(2,3) = THSO#(-.5*#R/RLAM) 
TT(2,4) = TH302®(-.5*R/RLAM) 
c 
TT(3,3) = PHI31/0SIN( ANG) 
TT(3,4) = PHI41/DSIN(ANG) 
TT(4,3) = TH&41®.5*R/RLAM 
TT(4,4) = TH31*.5*R/RLAM 
c 
15 TT(3,1) = PHI11/DSIN( ANG) 
TT(3,2) = PHI21/DSIN(ANG) 
TT(4,1) =-TH11*.5#R/RLAM 
TT(4,2) = TH21*.S2R/RLAM 
IF(IT.€O0.5) CALL ROWEX(TT) 
20 CALL TTINV(TT,HO) 
c 
IF(IFLAG.NE.O) RETURN 
IF(HO.EQ.0.) GO TO 30 
c 
TA(1,1) = TH20%(R*RLAM*#DSIN(ANGO)/(E*T)) 
TA(1,2) = THIO#(R*RLAMEOSIN(ANGO)/(E*T) ) 
TA(1,3) =-TH3O#(R*RLAM#OSIN(ANGO)/(E*T) } 
TA(1,4) = TH&4O=(R*RLAMZDSIN(ANGO)/(E*#T) ) 
TA(2,1) =-PHI20#(2.*RLAM*#2/(E*T)) 
TA(2,2) = PHI10*(2.*RLAM**2/(E*#T) ) 
TA(2,3) = PHI40*(2.*RLAM*2#2/(E*T)) 
TA(2,4) =-PHI302(2.*#RLAM#*#2/(E#T)) 
TA(3,3) =-TH31*(R*RLAM*DSIN(ANG)/(E 
TA(3,4) = TH41*(R*RLAM*DSIN(ANG)/(E 
TA(4,3) = PHI81*(2.*RLAM®22/(E*T)) 
TA(4,4) =-PHI314(2.*#RLAM**#2/(E*T) ) 

30 TA(3,1) = THZ21*(R*RLAM*DSIN( ANG) /(E*T)) 
TA(3,2) = THI1#(R*RLAM*OSIN( ANG) /(E*T) ) 
TA(4,1) =-PHI212(2.*RLAM**2/(E*T)) 
TA(4,2) = PHI11*(2.*RLAM**2/(E*T)) 
IF(IT.EQ.5) CALL ROWEX(TA) 

c 
DO 100 1=1,4 
BDO 100 J=1,4 
c=0.0 
oO 80 K=i,4 
80 CSC+TA(I,K)*TT(K,J) 
100 F(T, ,J)ec 
c 
RETURN 


END 


Cease Kee ARTETA RTT TAA ETA RTE T ERS ETERS EKER TESTES TESTE REESE ES 


SUBROUTINE CONE(IT,T,R,H,HO,E,PR,UW,F,TT,ANG,XM,RLAM, IFLAG) 


c THIS SUBROUTINE CALCULATES THE FLEXIBILITY MATRIX (F) FOR 
c & COMPLETE OR TRUNCATED CONE. THE FLEXIBILITY MATRIX IS 
c REDUCED TO A TWO BY TWO MATRIX FOR & COMPLETE CONE. 
c 
IMPLICIT REAL?8(4-H,0-2) 
DIMENSION F(6,6',TA(4,4),TT(4,4) 
(s 
ANG = R/57.295779513 
xM = (12.21. °PR¥*2))280.25 
RLAM= (XM*®24/(T*®DTAN( ANG) )**2)**0.25 
c 
x1 = 2.*RLAM*(H*®*0.5) 
KRIOKA SGT.91 19) 9) Go To 93939 
CALL MMKEL2(X1,BER21,BEI21,XKER21,XKEI21, 
*DBER21,0BE1I21,OKER21,0KEI21) 
FF "CHO.ONE O-.-) GO TO 10 
c 
c INITIALIZE THE MATRICES TO ZERO. 
oo 5 121,4 
po Ss J=i,4 
TANT, did) = “Oe 
Ss TT(I,JS) = O. 
Go To 15 
c 
10 KO = 2.*RLAM*(HO®#0.5) 
RECKOmar. 119.3 GO To 3398 : 


CALL MMKEL2(X0,BER20,BE120,XKER20, XKEI20, 


*DBER2 
tert AS 
Tories 
ere in hs 


nant 


©,DBEI20,0KER20,D0KEI20; 

Fy BER20/(-HO*#DSIN(ANG)?) 
BEI20//‘-HO*DSIN( ANG) ) 
XKER20/'-HO*OSIN( ANG) ) 


3) 
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SIN(ANG) ) 


®PR*BEI20)*T/(2.*KM2222HO) 
*PR*BER2Z0)*T/(2.*xXM**22HO) 
SPREXKEI20)*T/(2.*®XMe*22HO) 
#PREXKER20)#T/(2.*XM*#*22HO) 


N(ANG)) 
N( ANG) ) 
*PREXKEI21)#T/(-2.*xXM**2"H) 
=PR*EXKER21)*T/(-2.2xXM2*25H) 


*PR*BEI21)*T/(-2.*xXM2*22H) 
*PR*BER2Z1)*T/(-2.2XM*#24H) 
(TT) 


2.=PREBER20)/(2.*E*T/DSIN( ANG) ) 
2.=PREBEI20)/(2.*E*T/DSIN( ANG) ) 
2.=PR*XKER20)/(2.8E*T/OSIN( ANG) 
2.PR*¥XKEI20)/(2.8E*T/OSIN( ANG) 
2/xXMes2) 

2/XM*22) 

*2/XM*22) 

22/XME*2) 

=2/KXM222) 

=2/XM**2) 
2.=PR*XKER21)/(2.2E*T/DSIN( ANG) ) 
2.ePR*XKEI21)/(2.*E*T/DSIN( ANG) ) 
2.2PR*BER21)/(2.*E*T/DSIN( ANG) ) 
2.*PR*BEI21)/(2.*E*T/DSIN( ANG) ) 
2/XM*#*2) 

2/kM=22) 

(TA) 


) 
) 


TT(1,8) = XKEI20/(-HO#D 
c 
TT(2,1) = (XO*#DBEI20+2. 
TT(2,2) =-(XO*#OB8ER2042. 
TT(2,3) = (XO*DKET20¢2. 
TT(2,4) 2-(KO*#DKER20¢2. 
c 
TT(3,3) = XKER21/(H*DSI 
TT(3,4) = KKEI21/(H#OSI 
TT(4,3) = (X1#OKEI21%2. 
TT(4,4) =-(KX1*#OKER2142. 
c 
15 TT(3,1) = BER21/(DSIN( ANG) #H) 
TT(3,2) = BEI21/(DSIN( ANG) #H) 
TT(4,1) = (KX1*DBEI2142. 
TT(4,2) =-(X1*O0BER2142. 
IF(IT.EQ.6) CALL ROWEX 
20 CALL TTINV(TT,HO) 
c 
IF(IFLAG.NE.O) RETURN 
IF(HO.EQ.0.) GO TO 30 
TA(1,1) = (XO*DBER20 - 
TA(1,2) = (XO#DBEI20 - 
TA(1,3) = (XO*#DKER20 - 
TA(1,4) = (XO*DKEI20 - 
TA(2,1) = BEI20/(-EeTes 
TA(2,2) =-BER20/(-EeT#* 
TA(2,3) = XKEI20/(-EsT* 
TA(2,4) 2-XKER20/(-E*T# 
TA(&,3) = XKEI21/(-E*T# 
TA(&4,&) =-XKER21/(-E*Ts 
TA(3,3) = (K1*OKER21 - 
TA(3,4) = (K1*DKEI21 - 
30 TA(3,1) = (X13D0BER21 - 
TA(3,2) = (KX1#O0BEI21 - 
TA(4,1) = BEI21/(-E*T#% 
TA(4,2) =-BER21/(-E*T#s 
TF CFT 20.6) CALL ROWEX 
[se 
S ASSEMBLE THE ELEMENT FLEXIBILITY MATRIX F. 
c 
DO 100 12=1,4 
DO 100 Jv=1,4 
c=0.0 
DO 8O K=1,4 
80 CzC+TA(I,K)*#TT(K, J) 
100 F(I, J) s¢€ 
RETURN 


$39 WRITE(6, 
1000 FORMAT ( ‘ 
* . 
STOP 
END 
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1000) 
PROGRAM STOPP 
CHECK INPUT A 


ED FOR CONE’,/, 
RGUMENT FOR KELVIN FUNCTIONS’ ) 


SUBROUTINE MMKEL2(X,BER2,BEI2,XKER2,XKEIZ, 
*DBER2,0BEI2,DKER2,0KEI2) 


annan 


THIS SUBROUTINE CALCULATES KELVIN FUNCTIONS AND ITS 
DERIVATIVES OF THE FIRST AND SECOND KINDO, UTILIZING 
IMSL ROUTINES MMKELO AND MMKEL1. 


IMPLICIT REAL*8(A-H,0-2) 


CALL MMK 
CALL MMK 
R2 = 2.% 
BER2 = - 
BET2 = - 
XKER2 = 
XKEI2 = 


DBER2 
DBEI2 
DKER2 
DKEI2 


uwnn 


RETURN 
END 


CREST ESTES KAAS AT AAT SEK ARTS TTT SKK TT SSE TALS ET AE SKRAETTT KETC AERSEEAE 


SUBROUTI 
c THIS SUBROU 
c ALL TYPES O 
c 


ELO (X,BERO,BEI 
EL1 (X,BER1,BEI 
20.5 


R2/X * (BER1-BE 
R2/X * (BE11+BE 
-R2/X = (XKER1I- 
“R2/X = (XKEI1+ 


- (BER1+BE!I1)/R2 
-(BEI1-BER1)/R2 
- (XKER1+XKEI1)/ 
- (XKEI1-XKER1)/ 


NE TTINV(A,HO? 
TINE INVERTS TH 
F SHELLS 


©,XKERO, XKEIO, IERO) 
1,XKER1,XKEI1,1ER1) 


T1) -BERO 
Ri) -BEIO 
XKEI1) -XKERO 
XKER1) -XKEIO 


- 2.*BERZ/X 
2.*BEI2/x 

2.*XKER2/X 

2.*xXKEI2/XxX 


R2 
R2 


E (TT) MATRIX FOR 


IMPLICIT REAL*¥8(A-H,0-2) 


CIMENSIO 
Tt 


c 

3 
uunon 
B-b 


x 
owen 


Go TO 25 
10 CONTINUE 
15s CONTINUE 


20 DO 8oO J= 
po 75 J= 


N A(4,4),B(4,4) 


1,4 
3 gk 
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25 


30 
35 


40 


4s 
50 


60 
65 


70° 
75 
te) 
85 
go 


ss 
100 


Ni = 1 

B(I,J) = 0. 

DO 65 K21,4 
IF(K.EQ.1) GO To 65 
DO 60 L=1,4 
IF(L.EQ.J) GO TO 60 


N2 = 1 
G2) =o. 

DO SO M=11,12 
IF(M.EQ.K.OR.M.EQ.I) GO°TO SO 
DO 45 NeJ1,u2 
IF(N.EQ.L.OR.N.EQ.J) GO TO 45 


00 35 MMsI1i,12 

IF (MM.EQ.M.OR.MM.EQ.K.OR.MM.EQ.I) GO TO 35 
DO 30 NN=J1,J52 
TF(NN.EQ.N.OR.NN.EQ.L.OR.NN.EQ.J) GO TO 30 
c3 = A(MM,NN) 

IF(HO.NE.O.) GO TO 40 

B(M,N) = A(MM,NN) 

B(MM,NN) = A(M,N) 

GO TO 40 

CONTINUE 

CONTINUE 


C2 = C2 + (-1.)8*(N241)8A(M,N)*CO3 
N2 = N2 + 1 
IF(HO.EQ.0..AWND.N2.EQ.3) Go TO 8§ 
IF(N2.EQ0.3) GO TO 55 

CONTINUE 

CONTINUE 


B(I,J) = B(I,J) + (1.388 (N141)FA(K,L)*C2 
Ni = N11 + 1 

IF(N1.EQ0.4.A8ND.1.E9.1) GO TO 70 
IF(N1.E0.4.AND.1.GT.1) GO TO 75 

CONTINUE 

CONTINUE 


DET = DET + (-1.)#8(I+UDFA(1, J) 8B(I, J) 
CONTINUE 
CONTINUE 


GO To 90 
OET = C2 


DO 100 11,4 
DO 95 J=1,4 

A(I,J) = (-1.)##(1T4#JI)*B(JU,1)/DET 
CONTINUE 

RETURN 

END 
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20 


SUBROUTINE ROWEX(A) 


THIS SUBROUTINE PERFORMS ROW INTERCHANGES FOR 
MATRICES (TA) AND AND (TT) WHEN THE SHELL IS 
INVERTED. 


IMPLICIT REAL*8(A-H,0-2) 
DIMENSION A(4,4) 


TEMP = C¥A(I,J) 
A(I,J) = C#a(1+2,J) 
A(I1+#2,J) = TEMP 
Cres 

CONTINUE 

RETURN 

END 


CeeeeeeeeeeaeESESHS RSET ETT KSTKEE SES ERSEHKKKKAE THEATRE SET ETERS EKECEATSE* TEST 


c THIS SUBROUTINE COMPUTES CYLINDER PARTICULAR SOLUTION OISPLACEMENTS 


10 


SUBROUTINE PCYLIN(T,R,H,HO,E,PR,UW,ALPHA,F,PSO,IP,PV,N,PSF,PBT) 


IMPLICIT REAL*8& (A-H,O0-2) 
OIMENSION F(6,6),PS0(4),PSF(20,6) 


Do 10 1=1,4 

PSO(1I)=0.0 

PRiGrec Ure 1eORS VP. Gt tae Go TO 3939 
PSD(1)=2-PBT*PR#R/(E*T) 
PSO(3)=PS0(1)? 


GOL TFOe tZo, 30,267,450, 70770, SO EP 
PSD(1)=PV#R¥*¥2/(E*T) +PSD(1) 
PSD(3)=PS0(1)! 

GO TO 70 


PSD(3)=PSD(3)-UW*T*PR*R*H/ (EXT) 
C=UWtTFPR#R/(E*T) 
PSD0(2)=PSD(2)-C 
PSD(4)=PSD(4)-C 

Go TO 70 


=-ALPHA#R*PV 
PSD(1)=PSD(1)+C 
PSO(3)=PSD(3)+C 


GO TO 70 

C = PV#R#*2/( EFT) 
PSD(3) = PSD(3) + CH 
PSO (2). 2aPS0(2) <2€ 
PSD(4) = PSD(4) + € 
pO 100 1=1,4 
C=PSD(1) 

po 8o J=1,4 
C=C+F(1,J)*PSFI(N,J) 
PSD(1)=C 

RETURN 


WRITE(6,1000) IP 
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CERES SHEET RRS SARE KEK TRESEKRERERE SEER ERRATA KEE RERRSEERTEREES 
SUBROUTINE PDOME(IT,T,R,H,HO,E,PR,UW,ALPHA,F,PSD,IP,PV,IDV,PSF, 


FORMAT (°* PROGRAM STOPPED FOR CYLINDER IP =’,14) 
STOP 
END 


s ANG,ANGO,WHT,PBT,IFLAG,N) 


C THIS SUBROUTINE COMPUTES DOME PARTICULAR DISPLACEMENTS (PSD) 


nann 


15 


20 


25 
30 


35 


7° 


100 


c 
999 
1000 


CLASES ERE TESA ESET TTS TAKES TAETA STE KSA TASER SKELTER E REE EERE EEL ETK TERETE 


15 


50 


IMPLICIT REAL*8(A-H,O0-2Z) 
DIMENSION F(6,6),PS0(6),PSF(20,6) 


IF(IP.LT.1.0R.IP.GT.7) GO TO $939 
DO 10 121,6 

PSD(1I)20.0 
RLAM=(3.4#(1.-PRe*2)2(R/T)**2) 22.25 


PHI1 = ANGO 

PHI2 = ANG 

Ci 2 4. 

eS Orr Eo. 2) Go TO 15 
PHI1 = ANG 

PHI2 = ANGO 

Gima <7, 


INITIALIZE STRESS RESULTANTS AT THE TOP AND 
BOTTOM OF THE SPHERICAL SEGMENT. 


IF(PHI1.E0.0.) GO TS 20 
RN10 = FNI(PV,PST,UW,R,T,IP,IT,PHI1,ANGO,ANG,WHT, IDV) 
RN20 = FN2(PV,PBT,UW,R,T,IP,IT,PHI1,ANGO,ANG,WHT, IDY)} 
IF (PHI2.E0.0.) GO TO 25 
RN11 = FNI(PY,PBT,UW,R,T,IP,IT,PHI2,ANGO,ANG,WHT, IDV) 
RN21 = FN2(PV,PBT,UW,R,T,IP,I1T,PHI2,ANGO,ANG,WHT, IDV) 


GOTO(30,30,30,40,70,30,30),IP 


IF(PHI1.E90.0.) GO TO 35 
PSD(1) = PSD(1) -R*DSIN(PHI1)*(RN2O-PR*RN10)/(E#T) 
IF(IP.EQO.2) PSD(2) = PSO(2)-UWER2=((1.%PR)*(C1*O0COS(PHI1)**2 
= -1.)/OSIN(PHI1)-C1#*OSIN(PHI1))/E 
IF(IP.EQ.6) PSD(2) = PSD(2) - PV#R*((1.+PR)*(C1*DCOS(PHI1)*#2 
* -1.)°28C1#O0SIN(PHI1)*#*2)/ 
= (DTAN(PHI1)*#E#T) 
IF(IP.E&Q.7) PSO(2)=SPSD(2) - ((RN10-RN20)4(1.+PR) 
x +PVER*((1.+4PR)=(WHTtC1#R) 
* +C1#*R*OSIN(PHI1) 
= +C1*#(1.+PR)*R*OCOS(PHI1) 
= ))/(E*T*OTAN(PHI1)) 
IF(PHI2.EQ.0.) GO TO 70 


PSD(3) = PSD(3) -R*OSIN(PHI2)2*(RN21-PR#RN11)/(E*T) 
IF(IP.EQO.2) PSD(4) = PSD(4)+*UW*R*((1,%PR)*(C1*O0COS(PHI2)*#2 


* -1.)/DSIN(PHI2)-C1*OSIN(PHI2))/E 

IF(IP.EQ.6) PSD(4) = PSD(4&) + PV*R*¥DCOS(PHI2)#((1.+PR)2(Ci*® 
= DCOS(PHI2)**2-1.)°2.*C1*OSIN 
= (PHI2)*#2)/(OSIN(PHI2)*E*T) 


IF(IP.EO.7) PSD(4)=2PSD(4) + ((RN11°RN21)4(1.+PR) 
-(C1*#PVER*((1.+PR)*(WHT4C1*R*OCOS 


=x 

= (PHI1))*(DSIN(PHI1)/DSIN(PHI2) )**2 
= > C1#*(1.+PR)FR*(OCOS(PHI2) 

= + 2.*(0COS(PH1I2)**3-DCOS(PHI1)2*3)/ 
= (3.*O0SIN(PHI2)*®*2))))) 

x /(E*T#OTAN(PHI2) ) 

= + PVFR*R*DSIN(PHI2)/(E*T) 

GO TO 7° 


PSD(3)=PSD(3)-ALPHA*PV#R*¥DSIN(PHI2) 
PSD(1)}=PSD(1)-ALPHA*®PV*R*OSIN(PHI1) 


po 100 1=1,4 
c=PSD0(1) 

DO 8O J=1,4 
C=aC+F(1,J)*PSFI(N,J) 
PsD(I)=C 

RETURN 


WRITE(6, 1000) TP 

FORMAT(‘’ PROGRAM STOPPED FOR DOME IP =’,14) 
RETURN 

END 


SUBROUTINE PCONE(IT,7T,R,H,HO,E,PR,UW,ALPHA,F,PSO,IP, 
* PV,10V,PSF,ANG,WHT,PBT,IFLAG,N? 


c THIS SUBROUTINE COMPUTES CONE PARTICULAR DISPLACEMENTS (PSD) 


IMPLICIT REAL#8(A-H,0-7Z! 
DIMENSION F(6,6),PSD(6),PSF(20,6) 


YRUUP SLT 2th OR Po Giness Go To 993 

DON TO! W416 

PSD(1)=0.0 

V2" =H 

Yi = HO 

Ci = Me 

RE CET NG. 6) Go TO 15 

Goreme eat 

Yrs: h 

Y2 = HO 

C = .S*PR*(H**24+HO*#2)/Y22*2 

GOTO (20,30,20,60,70,40,50),IP 

DERO = PV#(PR-DTAN(ANG) ) 

DER1 = PV*(C -DTAN(ANG) ) 

Go TO 56 

DERO = UW*T*(PR/DCOS(ANG)-C1*D0SIN( ANG) *DTAN( ANG) ) 
DER1 = UWtT2(C /DCOS(ANG)-C1*DSIN( ANG) *DOTAN( ANG) } 
Go TO S6& 

DERO = PV*DTAN(ANG)?(PR-C1*DSIN(ANG)**#2) 
DER1 = PV*OTAN(ANG)*(C -C1*DSIN( ANG) **2) 
GOreT.0" 36 

DERO = PV*DTAN(ANG)*(-WHT-C1*2.*Y1*DCOS(ANG) 
2 + 0.5*PR*WHT#(1.+4+(HO/Y1)**2) 
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* + C1*2.2#Y1*0COS(ANG)/3. 

= + Ci®PR*HO**320COS (ANG) /(3.2#Y122#2) ) 
Si DERO = PV*DOTAN(ANG) *(-WHT-C1*#2.*Y22DCOS( ANG) 

* + O.S*PR*EWHT#(1.¢(HO/Y2)**2) 

z + C1*#2.*V2*0COS(ANG)/3. 

* + C1*PR*HO#*3*DCOS (ANG) /(3.*Y2*22) ) 
S56 IF(Y¥1.E9.0.) GO TO 57 


RN10 = FN3(IP,IT,Y1,ANG,H,HO,WHT,PV,PBT,UW,T, IDV) 
RN20 = FN4(IP,IT,Y¥1,ANG,H,HO,WHT,PV,PBT,UW,T) 
PSO(1) = PSO(1) - Y1*#OSIN(ANG)*(RN20-PR#*RN10)/(E#T) 
PSD(2) = PSO(2)-OTAN( ANG) *((1.+PR)*(RNI1O-RN20) 

* -~Y1*DERO)/(E#T) 

S57 PBB = 0. 

RN11 © FN3(IP,IT,Y2,ANG,H,HO,WHT,PV,PBB,UW,T, IDV) 
RN21 = FN4(IP,IT,Y2,ANG,H,HO,WHT,PV, PBB, UW,T) 
PSO(3) = PSD(3) - Y2*OSIN(ANG)*(RN21-PR#ERN11)/(E#T) 
PSD(4) =. PSD(4)-DTAN(ANG)#((1.+PR)*(RN11-RN21) 

= i -Y2*OER1)/(E#T) 

Go TO 70 


60 PSD(3)2PS0(3)-ALPHA#PV*Y2*0SIN( ANG) 
PSD(1)=PSD0(1)-ALPHA#PV*Y1*®0SIN( ANG) 


c 
70 DO 100 121,48 
C = PSD(I) 
DO 80 J=1,4 
80 C=C+F(1,J)*#PSF(N,J) 
100 PSD(I)=C 
RETURN 
c 


sss WRITE(6, 1000) IP 
1000 FORMAT(’ PROGRAM STOPPED FOR CONE IP 2’,14) 
RETURN 
END 
CHESS T ASSET SASS ESATA ERSTE KERR TAKES ASARATRER ERASER TKARETREREER EXE 
SUBROUTINE SOL(4,8,NN,NEQ) 
THIS SUBROUTINE SOLVES A SET OF LINEAR ALGEBRAIC EQUATIONS 
OF THE FORM A*X=B BY GAUSSIAN ELIMINATION, WHERE ‘A’ IS A 
SQUARE MATRIX. B IS THE RIGHT-HANOD SIDE VECTOR 
ON ENTRY. BUT IS OVERWRITTEN WITH THE SOLUTION VECTOR ’X’ 
DURING BACK SUBSTITUTION : 
IMPLICIT REAL*8(4-H,0-72) 
DIMENSION A(NN,NN),B(NN) 


nannn 


NL=NEQ-1 

BDO 250 N=1,NL 
IF(A(N,N).LE.O.) GO TO 500 
NISN+1 


00 100 J=N1,NEO 

100 A(N, J)SA(N,J)/AIN,N) 
B(N)=B(N)/A(N,N) 
BDO 250 JI2N1,NEQ 
IF(A(1,N}.EQ.0.) GO TO 250 
C=A(1I,N) 
DO 200 J=N1,NEQ 

200 A(I,J)2A(1,J5)°CFA(N, J) 
B(1)sB(1)-C#B(N) 

250 CONTINUE 

c BACK SUBSTITUTION 
MENEO 
B(M)=B(M)}/A(M,M) 
oc 400 N=i,NL 
M1=M 
M=M-1 
00 400 J=mM1,NEQ 

400 B(M)=B(M)-B8(J)*A(M, J) 
GO TO 600 

soo WRITE (6,1000) N 


CALL EXIT 

1000 FORMAT(’ ZERO OR NEGATIVE ELEMENT ON MAIN DIAGONAL OF TRIANGULARIZ 
1ED STIFFNESS MATRIX ¢ Ug FOR EQUATION NUMBER ‘’,14) 

600 RETURN 
END 


CES KSEE EKER ESKER EE ERE SAK SSE TERR TESTER EETRAEACAR ERE ERTIES 


SUBROUTINE BASE(IFLAG,T,R,H,HO,E,PR,UW,BB,S,TT,D) 


C THIS SUBROUTINE COMPUTES THE FLEXIBILITY MATRIX (Si 
C FOR A BASE SEGMENT ON AN ELASTIC FOUNDATION, OR 
C (IF IFLAG=1) THE MATRIX BB TO DETERMINE INTERNAL 
C DISPLACEMENTS AND STRESS RESULTANTS 
IMPLICIT REAL*8(A-H,0-2) 
DIMENSION S$(6,6),TT(4,4),B(4,4%),G6(4,4) ,BB(4,4) 
DIMENSION PHI(4),PHIP(4),PHIDP(4),PHITP(4),I1VEC(4) 
DATA IVEC/2,5,4,6/ 
c 


c FUNCTION DEFINITIONS 
FII(RO,RIJ=CM*RI/RO®*2+CP/RI 
FOI(RO,RI!=2.*C/RO 
FOO(RO,RI)=CM*RO/RI*¥*2+CP/RO 
FIO(RO,RI)=2.*C/R1I 


LM=2 

IF(HO.EO.0.) LM=1 

LN=2*LM 
D=E*T44#3/(12.=(1.-PR**2)) 
STIFL=(D/R)*¥*0.25 
IF(IFLAG.EQ.2) GOTO 300 


SIGN=1.0 

RD=H 

00 50 J=1,6 

Do S50 1=1,6 
so $(1,J)#0.0 

po §1 =1,4 

DO 51 J=1,4 
51 B(I,J)=0.0 


DO 100 I=1,LM 

T1=22(J-1)4%1 

I12=1141 

IF(1I.EQ.2) RO=HO 

PFLISEOS 2) SIGN=-1.0 

RDI=1./RD 

RDI2=RDI**2 

CALL BSHAPE(STIFL,RO,PHI,PHIP,PHIOP,PHITP? 
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c 

c abdo 


320 


c 
c 


CESSES SRTTER ARTETA TRS SETARERERA THERA TAK TARAKEKKTEKA SKA EKARKESERSTTKEKRE 


Cees 
c FOR 


c FOR 


c FOR 


nn 


FOR 


an 


FOR 


= 


* 


= 


= 


c 


CESSES KSEE EEK TESTERS SESE TAKE AERA ASAT SETTER EEE KAK KEEL TEKST EERE 


c THI 
c AND 


DO 80 J=1,LN 
B(12,45)=-02(PHITP(J)+RDI*PHIOP(J)-RDI2Z*PHIP(J))#*SIGN 
B(11,J3)20*(PHIDP(J)+PR*ROI*SPHIP(J))*SIGN 

G(12,J)2PHI (J) 

G(11,J)SPHIP( J) 

CONTINUE 

CONTINUE 


CALL JINVER(B,TT,4,LN) 
IF(IFLAG.EQ.1) GOTO 210 


00 200 I=51,LN 
ITSIveéc(!) 

DO 200 J=i,LN 
JJZIVEC (J) 

c=0.0 

DO 150 K=1,LN 
C=C+G(1,K)*TT(K, J) 
S(II,JjJJd)sC 


FLEXIBILITIES FOR IN PLANE STIFFNESSES 
IF(LM.EQ.1) GOTO 2¢5 
C=(H*HO)2#2/( Tt (H**2°HO#=2)*E) 
CM=(1.-PR)#C 
CP=(1.+PR)#C 
$(1,1)2F00(H,HO) 
$(1,3)=2FOI(H,HO) 
$(3,1)SFIO0(H,HO? 
$(3,3)=FII(H,HO) 

GOTO 210 
$(1,1)8H2(1.-PR)/( ET) 


RETURN 


CALL BSHAPE(STIFL,H,PHI,PHIP,PHIDP,PHITP) 
DG 320 J=1,LN 
BB(1,J)2D0*(PR*PHIDP(J)*PHIP(J)/H) 
BB(2,J)=0* (PHIDP(J)+PR#PHIP(J)/H) 
BB(3,J)2-D*(PHITP(J)+PHIDP(J)/H-PHIP(J)/H**2) 
BB(4,J)=PHI (J) 

XSHO/STIFL 

VPCRALT 2. Ke2.6 

HO=X*STIFL 

CONTINUE 

GOTO 210 

END 


SUBROUTINE BSHAPE(STIFL,RO,PHI,PHIP,PHIDP, PHITP) 


S SUBROUTINE EVALUATES THE PHI VECTOR AND ITS DERIVATIVES 


A BASE ON ELASTIC FOUNDATION SEGMENT 
IMPLICIT REAL*8(A-H,0-2Z) 

DATA RT,PI/2.0,3.1415926536/ 

DIMENSION PHI(4),PHIP(4),PHIDP(4),PHITP(4) 


P8=PI/8. 
CD1=1./(O0SORT(RT)*STIFL) 
SIG=RDO*CO!1 

RSIG=SDSORT(SIG) 
COSP=DCOS(SIG+P8) 
SINPSOSIN(SIG+P8) 
COsSmM=OCOS(SIG-P8) 
SINMSOSIN(SIG-P8) 
ETASRT#*O.75*DSORT(PI) 
CPHIP=OEXP(S1G)/(ETA#RSIG) 
CPHIM=PI*DEXP(-SIG)/(ETA*RSIG) 
CoO2z=CDi**2 

cbo3=CO2*CD1 


M PHI VECTOR 

PHI(1)=2CPHIP*COSM 
PHI(2)=CPHIP*SINM 
PHI(3)=CPHIM*COSP 
PHI(4)=CPHIM=SINP 


M PHIP VECTOR 

SIG2}I=1./(2.#S1G) 

SIGP=i1.+SI1G21 

SIGM=1.-S1IG2I 
PHIP(1)2CD1*(SIGM*PHI(1)-PHI(2)) 
PHIP(2)=CD1*(SIGM*PHI(2)+PHI(1)) 
PHIP(3)=-CO1*(SIGP*PHI(3)+PHI‘4)) 


PHIP(4)=CD1#(-SIGP*PHI(4)+PHI(3)! 

M PHIDP VECTOR 

C2=1./(2.*S1GF2#2) 
PHIDP(1)2=CO2*C2*PHI(1)+(SITGM*PHIP(1)-PHIP(2))*#CO1 
PHIDP(2)=CD2*C2*PHI(2)+(SIGM=PHIP(2)+PHIP(1))2CDI 
PHIDP(3)=CD2*C2*PH1(3)-(SIGP*PHIP(3)+PHIP(4))#CD1 
PHIDP(4)=CDO2*C2*PH1I(4)-(SIGP2PHIP(4)-PHIP(3))*CD1 
M PHITP VECTOR 

C2=2.*C2 

C3=-1./(SI1G#¥*3) 


PHITP(1)=CD3*C32*PHI"(1)+CO2*C2*PHIP(1)+(SIGM*PHIDP(1)- 


PHIOP(2))*CO1 


PHITP(2)=2CDO3*C3*PHI(2)+CO2*C2*PHIP(2)+(SIGM*PHIOP(2)+ 


PHIOP(1))#CD1 


PHITP(3)=CD3*C3*PHI(3)+CD2*C2*PHIP(3)-(SIGP=PHIDP(3)+ 


PHIDP(4))*CDi 


PHITP(4)=CD3*C3*PHI(4)+CD2*C2*PHIP(4)-(SIGP*PHIDP(4)- 


PHIDP(3))*CD} 


RETURN 
END 


SUBROUTINE JINVER(A,B,NDOIM,NEOQ! 


S SUBROUTINE INVERTS THE MATRIX A BY THE JACOBI METHOD 


STORES THE RESULT IN E 
IMPLICIT REAL#8(A-H,OQ-2' 
OIMENSION A(NOIM,1),B(NOIM,1) 
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c 
c INITIALIZE THE B MATRIX 
BO 100 J=1,NDIM 
00 100 I=1,NDIM 
100 B(I,J)s0.0 
DO 110 J=1,NEQ 
110 B(J,J)#1.0 
c 
C BEGIN JACOBI REDUCTION OF MATRIX A AND ALSO OPERATE ON B 
DO 600 N=1,NEO 
IF (DABS(A(N,N)).LT.1.00-06) GO TO 998 
C=1./A(N,N) 
Ni=SN¢1 
IF(N.EQ.NEQ) GO TO 410 
00 400 J=N1,NEQ 
AJ=A(N,J)8C 
BJSB(N,J)#C 
00 300 121,NEQ 
A(1,J)EA(I, J) -AJ#FA(I,N) 
300 B(I,J)=B(1,J)-BU*A(I,N) 
A(N,J)2AaJ 
400 B(N,J)=BJ 


410 DS SOO J=1,N 
BJ=B(N,J)*C 
BDO 450 121,NEQ 
450 B(I,J)#B(1,d0)-BUFA(1I,N) 
500 B(N,J)=Bu 


Cc 

600 CONTINUE 
RETURN 

c 

ss9 WRITE(6, 1000) N 

1000 FORMAT(‘O’,’ ** ZERO ELEMENT ON MAIN DIAGONAL FOR EQUATION’, 

* 14° INDICATES MATRIX IS SINGULAR’) 

stop 
ENO 

c 

c 


Cee eee eee ees KSAT ATE RKATAKEK SERRATE ETHERS ERASE TE TEER ES 
SUBROUTINE PBASE(T,R,H,HO,E,PR,ALPHA, UW,F,PSO,IP,PV,N,PSF, PBF ) 
c THIS SUBROUTINE EVALUATES THE PARTICULAR SOLUTION 
c DISPLACEMENTS FOR A BASE SEGMENT 
IMPLICIT REAL*8(A-H,0-2) 
DIMENSION F(6,6),PS0(6),PSF(20,6),PBF(6,20) 


c 
B00 10 121,6 
10 PSD(1I)=0.0 
S 
€ SELECT LOAD TYPE 
PRGrip oii. 1 ORs LP). Gti) Go Tro 3s¢e 
GO TG (20,30,40,50,70,70,60) ,I1P 
c 
c INTERNAL PRESSURE 
20 PSD(5)=PV/R 
PSD(6)=PSD(5) 
GO TO 70 
= 
c DEAD LOAD 
30 PSD(5)sUWsT/R 
PSO(6)=sPSD(5) 
Go TO 70 
c 
c IN-PLANE PRESTRESS 
40 PSD(1)=PV#H/T 
PS0(3)=PSD(1)*HO/H 
Go TO 70 
( 
c UNIFORM THERMAL 
Xe) PSD(1)=-PV*ALPHA#H 
PSD0(3)=-PV*ALPHA*#HO 
GO TO 70 
cS 
c LIQUID PRESSURE 
60 PSD(S)=PV#H/R 
PSD(6)=PSD(S) 
& 
7° DO 100 1=1,6 
“  c=Psd(I) 
DO 8O J=1,56 
8° C=C+F(1,J)*PSFI(N,J) 
100 PSD(I)=C 
RETURN 
c 
$99 WRITE(6,1000) 
1000 FORMAT(’O’ ,‘ *= PROGRAM STOPPED IN SUBROUTINE PBASE FOR DIAGNOSE’ 
x ,‘'D ERROR’) 
c 
STOP 
END 
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